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1.  INTRODUCTION 


The  Naval  Underwater  Systems  Center  (NUSC)  sponsored  Independent 
Research,  "Finite-Difference  Solutions  to  Acoustic  Wave  Propagations,"  has 
been  successful.  As  a  result  to  date,  a  useful  product  —  the  Implicit 
finite  difference  ( I FD)  computer  software  program  for  the  solution  of  para¬ 
bolic  equations  —  has  been  developed  for  research  and  application  purposes. 
This  software  Is  now  being  used  Internationally  In  a  number  of  research 
laboratories  as  well  as  universities.  In  relation  to  the  development  of  the 
software  package,  the  theoretical  development  attracted  a  number  of  inter¬ 
nationally  well-known  scientists.  In  1982,  the  Office  of  Naval  Research 
(ONR)  Mathematics  Group,  under  the  coordination  of  Dr.  Richard  L.  Lau, 
awarded  a  research  grant  to  NUSC  to  encourage  technical  collaboration  with 
university  scientists  at  the  Yale  University  Center  for  Scientific  Computa¬ 
tion.  These  developments  set  the  stage  for  four  visiting  scholars  to  spend 
the  summer  of  1983  performing  research  aimed  at  the  solution  of  underwater 
acoustic  wave  propagation  problems  In  all  dimensions  (mathematically, 
physically,  and  computationally). 


This  report  is  arranged  In  sections.  Each  section  reports  the  tech¬ 
nical  accomplishments  for  a  particular  combination  of  authors.  Some  compu¬ 
tations  were  performed  by  VAX  11/780  computers  both  at  NUSC  and  at  the  Yale 
University  Computer  Science  Department. 


The  four  visiting  scholars,  all  professors  at  their  respective 
Institutions,  were 


Frederick  0.  Tappert 
University  of  Miami 
Miami,  FL  33149 

William  L.  Selgmann 
Rensselaer  Polytechnic  Institute 
Troy,  NY  12181 


Gregory  A.  Kreigsmann 
Northwestern  University 
Evanston,  IL  60201 

Donald  F.  St.  Mary 
University  of  Massachusetts 
Amherst,  MA  01003 


Other  academic  contributors  to  the  work  reported  here  were 


Martin  H.  Schultz 
Yale  University 
New  Haven,  CT  06520 


Kenneth  Jackson 

University  of  Toronto 

Toronto,  Ontario,  Canada,  MSS  1A7 


Navy  contributors  to  this  document  were 

Henry  Weinberg  Ding  Lee 

NUSC,  New  London,  CT  06320  NUSC,  New  London,  CT  06320 
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2.  A  WIDE  ANGLE  THREE-DIMENSIONAL  PARABOLIC  WAVE  EQUATION 


William  L.  Slegmann 
Rensselaer  Polytechnic  Institute 

Ding  Lee 

Naval  Underwater  Systems  Center 

Gregory  A.  Krlegsmann 
Northwestern  University 


ABSTRACT:  A  simple  extension  of  the  standard  two-dimensional  para¬ 
bolic  wave  equation  to  the  three-dimensional  case  can  be  accomplish¬ 
ed  by  retaining  the  angular  derivative  term.  This  extension  Is 
limited  to  dealing  with  small  vertical  angles  of  propagation.  A  new 
wide  angle,  three-dimensional  partial  differential  equation  Is 
developed  to  predict  the  sound  propagation  In  a  three-dimensional 
ocean.  This  formulation  Is  achieved  by  operator  theory  whose  mathe¬ 
matical  derivation  is  given  In  detail.  The  validity  of  the  formula¬ 
tion  Is  examined  In  full  through  discussion  of  approximation  and 
multiple  scale  analysis. 


2-1 /2-2 
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INTRODUCTION 


A  three-oimensional  parabolic  equation  (PE)  was  aevelopeci  by  Tappert^- 
almost  a  uecaae  ago.  It  was  the  same  equation  Tappert  usea  to  derive  the 
two-dimensional  paraoohc  wave  equation,  recognized  as  the  stanoaro  small 
angle  PE.  Encountering  the  three-dimensional  effect  in  the  ocean,  Baer^ 
initiatea  the  application  of  the  primary  three-dimensional  PE  to  real 
pruulems.  Recently  Perkins  and  Baer  implemented  the  Split-step  algorithm0 
into  a  computer  coae  to  solve  the  three-dimensional  PE.  Application  of  this 
tnree-dinierisioriai  coae  demonstrated  success  in  solving  three-dimensional 
proolems.  Prior  to  Baer  anu  Perkins'  three-dimensional  applications, 
pierce^  formulated  a  simplified  three-dimensional  parabolic  wave  equation 
expressing  one  spatial  variaole  in  terms  of  arc  length.  It  is  seen  from  the 
extension  of  tne  two-oimensionai  stanoara  PE,  the  Tappert  three-uimensiona I 
PE,  implemented  by  Baer-Perkins  for  real  applications  (for  simplicity  we  refer 
to  tne  equation  as  the  3D  PE),  only  handles  the  small  angle  propagation.  So, 
Pierce  has  not  pursued  nis  development  further.  It  is  the  purpose  of  this 
paper  to  report  tne  development  of  the  wide-angle  three-dimensional  PE,  wnicn 
accommodates  the  3D  PE.  During  the  course  of  this  wide  angle  development,  a 
riumuer  of  practical  questions  arose.  We  hign light  the  importance  of  these 
questions  ana  try  to  answer  these  questions  reasonably.  The  motivation  of 
answering  these  questions  led  us  to  the  formulation  of  the  three-uimensiona I 
wide  angle  PE.  These  questions  help  to  define  the  region  of  validity  and 
suggest  when  anu  where  the  three-dimensional  problem  can  ue  solveu 
two-dimensional ly.  A  formulation  based  on  the  operator  theory  is  a  starting 
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point;  the  complete  detail  is  cnscusseo.  An  analysis  using  the  multiple  scale 
technique  is  included  to  justify  the  operator  formalism.  A  selected  exact 
solution  has  been  usea  by  Schultz  et  al.5  to  discuss  the  validity  of  the 
formulation  as  well  as  the  accuracy  of  the  solution.  In  this  paper,  a 
simulated  tnree-oiniensional  problem  ana  an  application  are  included  to 
demonstrate  the  three-dimensional  wide  angle  PE  capability.  All  computations 
were  performed  on  the  VAX  11/780  computer  using  the  Vale  Sparse  technique.^ 


OPERATOR  FORMALISM 


We  oegir;  from  the  three-dimensional  Helmholtz  equation  for  the 
spatially  varying  part  of  the  acoustic  pressure  p  =  p(r,e,z),,  written  here  in 
cylindrical  coordinates  (r,o,z),  i.e., 


The  complex  pressure  is  p  times  e-i^t,  where  w  is  acoustic  frequency  in 

raa/s.  In  Eq.  (1),  kQ  _  Ui/Cq  ^na  the  index  of  refraction  is  n  =  n(r,e,z) 

=  co/c,  where  c  =  c(r,»,z),  the  oceanic  sound  speeu,  ano  Co  is  a  reference 
sound  speed.  A  thorough  discussion  of  conditions  and  assumptions  under  which 

Eq.  (1)  applies  to  oceanic  sounu  propagation  has  been  given  by  Pierce. ^ 
Boundary  conditions  for  Eq.  (1)  are  to  be  specified  at  the  ocean  surface  ano 
bottom.  The  source  term  is  omitted  from  the  right  siue  of  Eq.  (1)  i ri 
anticipation  of  PE  approximations  that  are  valid  away  from  the  source,  which 


is  assumeu  near  r  =  0. 
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Following  Tappert,!  we  let 


p(r,o,zj  =  u(r ,o,z)  v(r) 


(2) 


in  whicn  the  factor  v(r)  represents  a  rapidly-varying  portion  of  the  pressure 

ana  u(r,o,z)  is  its  modulation.  Substituting  Eq.  (2)  into  Eq.  (1)  yielus 

2  2 
a  i 

az‘ 


d 

.ar 


t"  +  (A  +  l  *1)  ®!t  +  +  iUi  +  k2  n2  u 


ra2v  1  av 

17 


0 


(3) 


It  follows  from  Eq.  (3)  that  it  ari  oscillatory  function  v  is  aetermirieo  as  a 
solution  of 


♦  i  HI  + 

r  ar 


v  =  0 


(4) 


then  the  u  satisfies 

L"  +  (A  +  1  H)  iH  +  L 

.  2  'r  v  ar'  dr  z 
ar  r 


4  *  4  * 

ao  az 


1)  u  -  0 


The  outgoing-wave  solution  of  Eq.  (4)  is 

V(r)  =  Hji)(k0r)  , 


(b) 


(b) 


where  H^l)  is  the  Hankel  function  of  zero-th  oraer  of  the  first  kina. 
Since  the  parabolic  approximation  is  aesireu  for  the  solution  at  large 
distances  from  the  source,  it  is  appropriate  to  apply  a  furfielu 
approximation,  which  is  expresseu  by  kor  >>  ] #  oefer  until  the  next 
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section  a  discussion  of  tne  validity  ana  tne  quantification  of  this 
assumption.  For  now  we  employ  it  to  approximate  Eq.  (6)  by  an  asymptotic 
expansion 


kor 


(7) 


Using  Eq.  (7)  in  Eq.  (5)  gives 


8  U  _j_  Oil,  3U  1  6  U  3  U  ,, 2  / 


so 


+  ^  +  k0^n 

az 


i) 


u  =  0 


(&) 


If  the  first  term  in  Eq.  (8)  is  neglected,  we  obtain  upon  rearrangement  a 
fundamental  30  PE,  whicn  is  Eq.  (7.7)  of  Ret.  1,  i.e., 


i  k  2 

ju  0.2,  .  ,  ,  l  a  u  .  l  3  u 

—  =  -m—  [n  (r,e,z)  -  lju  +  - -  +  - 7  — 7 

3r  1  ^Ko  sz^  2kar^ 


(9) 


Neglecting  the  term — aria  regarding  n(r,u,z)  as  azimutnaliy  lriuepenuent 
2kQr  3<a 

the  stanaara  two-uiniensional  PE  results,  i.e., 


ik 

3  U  0 

a?  =  t~ 


[n2(r,z)  -  lju  + 


(10) 


Equation  (9)  nas  been  exploited  in  calculation  of  sounu  propagation  through  a 
•mesoscale  eaay.2  if  tne  last  term  in  Eq.  (9)  is  neglected  but  azimuthal 
dependence  is  retained  m  n(r,e,z),  then  a  simpler  PE  is  obtained  for  which  an 
efficient  implementation  has  been  demonstrated. 3  This  equation  is  useful 
specifically  in  the  absence  of  horizontal  diffraction  of  acoustic  energy,  as 
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for  example  with  weak  azimuthal  sound-speea  variations  ana  without  azimuthal 
redirection  of  energy  from  boundary  interactions.  Finally,  it  all  azimuthal 
dependence  is  neglectea  in  Eq.  (9),  the  usual  two-dimensional  PE  remains.  As 
is  well  known,  Eq.  (9)  arid  its  simplifications  are  valid  for  narrow  vertical 
angles  of  propagation,  in  order  to  obtain  a  3D  PE  appropriate  for  wider 
angles,  we  first  employ  an  operator  formalism. 


We  return  to  Eq.  (8)  ana  express  it  in  operator  form  as 

Pt  +  2lko  Ir  +  ^*7  +  T  ^7  +  ko^2  -  x)l  u  =  0 
(_a  r  3z  r  as> 


(11) 


An  approximation  to  Eq.  (11)  is  made  by  factoring  the  operator  as  follows: 

3  i  i  i.  J  I.  nl  ,  3 


lTF  *  -  ,Ko  «  Ljr  -  .k0  .  ik0  q:  U  .  0  , 


where 


Q  = 


l  +  ("2  -  1)  +  h1-! +-T71~I 

ko  3 7  kor  30 


1/2 


(12:) 


(13) 


Equations  (11)  and  (12)  are  not  equivalent  because  the  operators  Q  ana  a/dr  do 
not  in  general  commute.  However,  provided  these  operators  are  in  some  sense 
nearly  commutative,  it  is  appropriate  to  regard  Eq.  (12)  as  a  factorization 
approximation  of  Eq.  (11).  We  make  this  approximation  ano  will  aiscuss  later 
its  validity.  The  solution  of  Eq.  (12)  consists  of  waves  incoming  ano 
outgoing  in  the  raoial  directions,  ana  we  neglect  the  incoming  wave  (the 
second  factor  in  Eq.  (12)),  which  is  usual  in  the  PE  method.  Therefore,  the 
envelope  u(r,»,z)  satisfies  the  formal  equation 
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au 

ar 


1kou 


=  uo  Qu 


(14 


De term i nation  of  u  requires  some  definition  of  the  operator  Q  in  Eq.  (14). 


We  specify  Q  by  first  expressing  it  as 

Q  =  L1  +  X  +  Yj1/2  ,  (1! 

ana 

X  =  (n2  -  1)  +  i-ll  ,  Y  =  l~.  iL-  .  (K 

k2  az2  (kQr)2  ao 

The  funoamental  3D  PE, Eq.  (9),  can  oe  obtainea  by  expanding  the  square  root  in 
Eq.  (4b}  in  a  Taylor  series  ana  retaining  only  tne  linear  terms  in  X  ana  Y. 
Rather  than  a  (linear)  polynomial  approximation  for  Q,  we  use  a  rational 
function  approximation,  i.fe., 

„  .  1  *  P1X  *  P2V 

■>  *  1  v  q .K  *"'CT  •  1 

where  >  p^,  ,  aria  q^  are  constants  to  be  chosen.  The 

interpretation  of  the  fraction  in  Eq.  (17)  is  premultiplication  of  the 

numerator  by  the  inverse  of  tne  oenomiriator .  Tnus  when  Eq.  (i7)  is  inserteu 
in  Eq.  (14),  the  equation  governing  u  becomes 

+  ikQu  =  ikQ[l  +  qLX  +  q^YJ-1  [1  +  PjX  +  p^Yju  (If 

or,  equivalently, 

L1  +  q  X  +  q-Yj|^  =  ik  L(p ,-q ,)X  +  (p  -q  )Yju  .  (li 
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We  note  that  wnen  q^  _  q2  =  q,  [q.  (18)  reduces  to  the  yarrow-angle 
3D  PE  of  Eq.  (9)  for  the  values  =  p2  =  1/2,  which  are  just  those  in  the 
linear  Taylor  series  for  Q.  For  the  two-dimensional  problem  Y  =  0, 
rational-function  approximations  have  been  discussed. ?  In  particular,  the 

choices  p^  _  2/^,  =  i/4  for  the  two-dimensional  case  P2  =  q2  =  0 

nave  been  suggesteu  by  Claeroout^  for  wiaer-angle  propagation.  These  values 
are  precisely  those  necessary  for  an  approximation  to  Q  in  Eq.  (17)  correct  to 
quadratic  terms  in  X.  The  analogous  result  for  the  three-oimensiona i  case  is 
touno  by  squaring  Eq.  (17)  anu  matching  coefficients  of  X,  Y,  X^,  Y^,  ana  XY. 
It  can  be  shown  that  the  resulting  five  equations  are  satisfied  by  the  four 

choices  P(  =  p2  =  3/4  ano  qi  =  q2  =  1/4.  Thus,  these  values  give  a 
rational-function  approximation  to  Q  correct  to  second  order  in  the  operators 
X  ana  Y.  We  use  them  iri  this  paper  to  specify  a  wiaer-angle  3D  PE  from  Eq. 


(19),  i.e., 


i  .  1/-2  ,  ,  ,  1  32  ,  1  1  1  3U 

1  +  j(h  -  1)  +  — 7  — ?  +  - 7  —7 

L  4  4k2  az2  4(kQr)2  ae2J  9r 


ik  r  , 


1 )  ♦  -i 

ko  3Z  (ko 


1  a2  ] 

W  U 


1  d 

Neglecting  the  terms  involving  - 7  — *■  ano  regarding  n(r,e,z)  as 

(k  r )2  as2 

azimutnally  independent,  the  twu-aimensional  wide  angle  PE  results  in  the 
sense  of  ClaerDout 
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i  k 

Q 
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[n  (r,zj  -  1) 


1 

~2 - 7 

ko  « 


]• 


UD 


Note  that  Eq.  (20)  is  a  third-order  partial  differential  equation  and  a 
discretizeo  version  has  been  analyzed  for  numerical  stability. 9  Other 
choices  for  tne  parameters  plt  P2j  qit  anG  Q2  have  also  been 

investigated. 10 


DISCUSSION  OF  APPROX IMAT10NS 

The  wide-angle  3D  PE  of  Eq.  (20)  was  derived  subject  to  a  number  of 
assumptions  anti  approximations.  The  principal  advantages  and  limitations 
conmon  to  all  PE  approximations  are  oiscusseo  in  Ref.  1  (see  also  Ref  11). 

For  applications,  we  are  particularly  interested  in  determination  of 
limitations  on  oceanic  ranges  where  Eq.  (20)  is  appropriate  and  where  tne 
azimutna l-oerivative  terms  in  Eq.  (20)  are  significant.  We  focus  here  on 
three  of  the  assumptions  used  in  the  preceding  section:  farfield, 
factorization,  aria  rational-function  approximations.  Our  discussion  leads  to 
suggestive,  rather  than  rigorous,  conditions  specifying  range  intervals  where 
Eq.  (20)  or  simplifications  of  it  should  be  employed.  These  conditions  are 
supported  by  arguments  in  this  section;  asymptotic  derivations  and  numerical 
resu  I  ts  wi I  I  follow. 

We  turn  first  to  the  farfield  approximation.  Some  indication  of  the 
range  beyono  which  this  approximation  applies  is  very  desirable  in  order  to 
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provide  an  estimate  for  the  minimum  range  of  applicability  for  the  3D  PE. 
Denote  by  v(r)  the  right  sice  of  Eq.  (7),  i.e.,  the  leaning  term  in  the 
asymptotic  expansion  for  v(r)  in  Eq.  (6)  for  large  kQr.  Suppose  we  choose 
to  regaro  7  as  an  acceptable  approximation  to  v  if  the  relative  difference  in 
their  moduli  is  less  than  some  tolerance  «,  i.e.,  if 


I  v  ( r )  | 


I  v  ( r )  |  _1  <  6 


(22) 


This  condition  focuses  only  on  differences  in  mooulus,  rather  than  iricluoing 
differences  in  phase,  which  are  of  less  interest  here.  Now  it  is  known^ 
that 

I  v(r)|  -  |  v ( r ) |  1 - 1 - -  ♦  0  I— zr1!,  k  r  -  -  .  (23) 

I  16(k0r)'  L(kor)4J|  0 


In  Eq.  (23),  the  terms  in  braces  ({  })  alternate  in  sign  and  have  the  property 
that  the  remainder,  after  retaining  any  number  of  terms,  is  no  bigger 
than  the  first  term  neglected.  From  Eq.  (23)  it  follows  that  the  7  is 
regaroea  as  acceptably  approximating  v  if 


k  r  > 


4 *Jl 


(24) 


Iri  terms  of  acoustic  frequency  f,  Eq.  ^4)  requires  range  r  to  satisfy 


r  i  rf  =  8»fts/4 


(2b) 


in  whicn  r|  1S  minimum  range  for  the  farfield  approximation  to  apply. 
For  example,  suppose  6  =  0.01,  corresponding  to  differences  between  v  and  7 


2-11 


TD  714b 


being  bounaeu  by  1%.  Then  as  f  increases  from  10  to  200  Hz,  rf  ciecreases 
from  about  60  m  to  about  3  m.  An  alternative  expression  for  this  case  is 
*V  >  2.5. 


Adoption  of  a  criterion  such  as  the  above  suggests  neglect  of  terms  of 
appropriately  small  estimated  size  in  the  governing  equations.  For  instance, 

the  approximation  7(r)  satisfies 


d2v  1  dS  ,  2  - 

^  ^  k°  V 


1  -  - =  0  , 

4<k  r)2J 


(26) 


rather  than  Eq.  (4)  satisfied  by  v(r).  If  Eqs.  (24)  ana  (25)  hold,  then  the 
last  term  in  Eq.  (26)  is  no  bigger  than  4a.  Thus,  the  approximation  of  v  by  v 
is  tantamount  to  neglecting  tnis  term,  wrnch  for  6  =  0.01  is  of  relative 
magnitude  no  bigger  than  4«.  This  behavior  is,  of  course,  typical  of  a 
regular  perturbation  for  wmch  neglect  of  a  term  of  some  small  size  produces 
an  error  of  comparable  size  in  the  solution.  It  follows  that  unless  any  term 
in  a  governing  equation  is  capable  of  producing  a  singular-perturbation 
effect,  it  is  apparently  consistent  to  ignore  the  term  if  its  relative  size  is 
no  bigger  than  about  4%.  An  immediate  application  of  this  criterion  is  in  the 
simplification  of  the  coefficient  of  3U/ar  in  Eq.  (5).  Using  Eq.  (7)  and  the 
result  tnat  tne  asymptotic  expansion  of  dv/or  is  the  derivative  of  tne 
expansion  for  v,  it  is  easily  shown  that  the  coefficient  of  au/ar  in  Eq.  (8) 
is  multiplied  by  (1  +  l/8(k  r}2).  However,  inequality  (24)  means  that 
this  factor  in  the  farfield  approximation  is  no  bigger  than  (1  +  2s.)  Thus, 
it  is  apparently  consistent  to  ignore  tnis  factor  in  the  farfield 
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approximation.  The  inequality  (24)  with  «  =  0.01  (for  instance)  ensures  that 
no  more  than  4%  error  is  committed  in  terms  of  botn  the  u  anu  v  equations 
in  the  farfieid  approximation  and  that  the  error  in  the  modulus  of  7  is  even 
smaller.  It  can  be  shown  by  using  another  asymptotic  expansion  that  the  error 
in  the  phase  of  7  is  actually  of  comparable  magnituae. 


As  employed  in  this  paper,  the  rational-function  approximation  to  the 
square-root  operator  Q  is  given  by  Eq.  (17)  with  p1  =  P2  =  3/4  ana  =  q2  =  1/4. 


For  convenience  this  is  rewritten  as 

Q 

where 

Z  -  X  +  Y  *  n2 


1  ♦  (3/4)2 

iv  (1/412  * 


1  +  + 

k0  az 


1 


0 


(27) 


(28) 


As  mentioned  previously,  a  primary  advantage  of  this  approximation  is  that  it 
is  correct  to  second  oraer  in  Z,  i.e., 

I  -1  3 

Q  =  (1  +  *Z)  d  + 

=  1  ♦  \l  -  jZ2  +  0(Z3)  .  (29) 

An  alternative  expression  of  this  fact  is  that  tne  only  other  condition  needed 
for  Eqs.  (11)  and  (12)  to  De  identical,  in  addition  to  commutativity  of  Q  and 
a/ar,  is 


Q2  -  1  =  2 


(30) 
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From  Eq.  (29)  it  can  be  shown  that  Eq.  (30a)  holds  to  terms  of  0(Z3).  Next, 
we  note  that 


k  az  (k  r)  a  o 
0  0 


7? 

0 


4 

a 

~~2~1 

ao 3z 


111 L.. -J1 


r  2  2i 

3^  +  1  4 

'4 1”2  - 

1)  +  i_  (n*  -  i; 

)  , 

(31) 

7  77_ 

ko 

_8Z 

4 

CD 

s. 

- 

ana  it  follows  that  the  effects  of  the  terms  in  Eq.  (31)  are  incluoeo  by  Eq. 
(27).  Thus,  even  though  Eq.  (20)  explicitly  contains  no  fourth-order  2  and  » 
derivatives  the  effects  of  fourth-order  derivatives  in  Eq.  (31)  are  in  some 
sense  incorporated  properly  into  Eq.  (20).  On  the  other  hand,  Eq.  (20)  does 
not  contain  the  effects  of  any  sixth-order  aerivatives,  such  as  those 
appearing  in  Z3,  or  similar  terms  like  (n2  -  l)3  or  (n2  -  l)2 
(kQ)-2  (a2/a2z).  The  comparison  of  terms  neglected  with  those 
retained  is  most  easily  accomplished  by  scaling  and  asymptotic  expansions  such 
as  those  iri  the  next  section.  The  purposes  of  the  limited  discussion  here  are 
to  indicate  which  types  of  terms  are  modeled  correctly  by  the  approximation 
Eq.  (27)  ana  to  proviae  a  basis  for  examining  the  factorization  approximation. 


The  factorization  approximation  is  exact  when  Q  from  Eq.  (13)  and  a/ar 
commute.  Since  this  is  not  true  in  general,  Eq.  (12)  can  be  exparioeu  to  yielu 
Eq.  (11)  with  the  additional  term 


ik  I  ^ — (J  -  Or— j  u  =  u 
oLarN  ^ar-1 


(32) 
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on  the  left  side.  To  appraise  the  neglect  of  this  unbounded  operator,  we 
compare  its  terms  witn  those  retained  in  the  equation  for  u.  We  use  the 
expression  Eq.  (27)  for  Q,  but  other  definitions  from  Eq.  (29)  or  other 
parameter  choices  in  Eq.  (i7)  coula  be  treated  similarly.  In  view  of 
Eq.  (29),  it  follows  from  Eq.  (32)  that 


*  S(1)U  +  !^u  , 


where 


(33) 


\~~l  -  z  —j 

c,  ^  i-ir  L  * 


ar- 


(34) 


ana 


.(2)  lku  a  72  72  a 

*  '  T  lIF  Z  "  Z  Jr- ] 


(35) 


Using  Eq.  (28)  we  fine  that  the  leading  term  U)  has  the  form 


(36) 


As  with  the  farfiela  approximation,  a  comparison  should  be  made  here  of  terms 
neglecteo  (the  largest  of  which  are  U)u)  with  those  retained  Lin  Eq. 

( 1 1 ) J .  It  follows  tnat  the  factorization  reiies  un  kon3n/3r  being  small 

compared  to  k2(n2  -  1),  and  (k  r 3 )-l  (a2u/a©2)  being  small 

compared  to  r-2(a2(j/a©2) .  Since  n  is  close  to  one,  these  conditions  are 

k"1^2  -  l)'1  (an/ar)  small  (37) 

ana 

(hor)~^  smdl  1  ’  (38) 
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The  condition  [Eq.  (37)]  of  sufficiently  small  range  variation  of  the 
index  of  refraction  is  anticipated  from  analysis  of  the  validity  of  the 
two-dimensional  PE^ .  Equation  (38)  is  related  to  the  justification  of  the 
farfield  approximation  since  inequality  (24)  can  be  written  as  ( k Qr ) ~ ^  £46. 
Thus,  for  6  sufficiently  small  and  Eq.  (37)  valid,  both  farfield  and 
factorization  approximations  are  satisifed.  Furthermore,  Eq.  (38)  and  Eq.  (37) 
can  be  quantified  by  recalling  the  argument  following  Eq.  (26).  It  is 
apparently  consistent  to  ignore  the  effect  of  the  second  term  on  the  right  of 
Eq.  (34)  if 

(kQr)“1<46  .  (39) 

Here  we  have  already  neglected  terms  in  governing  equations  of  this  relative 
magnitude.  Equivalently,  for  (4  6  )-1  =  25  (when  6  =  0.01),  the 

second  term  in  Eq.  (36)  must  be  ignored;  for  f  *  10  Hz  (or  200  Hz),  this 
corresponds  to  ranges  bigger  than  about  600  m  (or  30  m) .  We  note  that  this 
represents  a  conservative  estimate  for  the  neglect  of  the  term,  which  may  in 
fact  have  an  insignificant  effect  for  even  smaller  ranges.  Also,  a  similar 
expansion  of  Eq.  (35)  and  a  comparison  of  terms  neglected  with  those  retained 
could  be  carried  out.  This  process  yields  Eqs.  (37)  and  (38)  along  with  other 
conditions  on  slowness  of  n(r,e,z)-variations  involving  various  partial 
derivatives  up  to  third  order  of  n(r,e,z).  We  omit  these  conditions  for 
brevity. 
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To  summarize,  the  farfiela  approximation  has  been  arguea  as  valid  for 

ranges  r  Digger  tnan  r^  given  by  Eg.  (2b).  Similarly,  the  factorization 
approximation  in  conjunction  with  the  rational-function  approximation  of  the 
operator  Q  is  appropriate  for  slow  variations  in  n(r,e,z)  |.see  Eq.  ( 37 ) j  ana 

for  r  bigger  tnan  «^2r^  L$ee  (39)j.  When  these  conditions  hold,  the 
three-dimensional  wide  angle  PE,  Eq.  (20),  should  be  applied.  As  range 

1  a^u 

increases  such  tnat  — = — 7  — 7  is  negligible,  the  two-dimensional  wide  angle 

(kQr)  so 

PE,  Eq.  (21),  should  be  used.  Iri  the  two-dimensional  application,  if  we 
regard  the  e-partitions  as  N,  this  coincides  exactly  with  the  "N  x  2D  Problem" 
defined  by  Perkins  ana  Baer. 3  it  is  important  to  remark  that  Eqs,  (37) 
through  (39)  are  not  predicated  on  any  statement  concerning  the  size  of  the 
o-vanation  of  u.  The  conditions  for  validity  of  both  the  farfiela  ana  the 
factorization  approximations  are  independent  of  whether  or  not  the 
o-oerivatives  in  Eq.  (20)  affect  the  propagation  significantly.  One 
resolution  of  this  latter  issue  is  provided  by  using  the  scaling  arguments 
presentea  below  to  compare  the  (kQr)_^  ( a2 / 302 )  terms  with 
32u /ar2.  It  is  sufficient  here  to  remark  that  azimuthal  variation  in  u 
can  be  introduced  by  three  mechanisms:  water-column  variations;  n  boundary 
fluctuations,  either  from  bottom  topography  ana  structure  or  from  surface 
irregularities;  and  directionality  in  the  representation  of  the  source 
nearf iela. 
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MULTISCALE  FORMALISM 


In  this  section  we  provide  a  systematic  asymptotic  derivation  for  a  class 
of  wide-angle  3D  PEs.  One  advantage  is  that  the  farfielo,  rational-function, 
and  factorization  approximations,  which  were  explicitly  required  in  the 
previous  oevelopment,  are  not  necessary  here.  Instead  the  scaling  ano 
asymptotic  expansions  produce  the  effects  of  these  approximations  without  any 
additional  assumptions.  Moreover,  considerable  insight  is  gained  into  the 
nature  of  these  approximations  and  the  conditions  for  their  validity.  The 
wice-angle  3D  PE,  Eq.  (30)  is  found  as  an  important  case  urioer  well-oefinea 
conditions. 


We  begin  with  the  Helmholtz  equation  LEq .  (1)J  which  does  not  incorporate 
source  effects,  in  this  paper,  we  do  not  treat  a  scaling  and  asymptotic 
expansion  appropriate  to  the  near-source  region.  Consequently,  we  regard  the 
spatial  portion  of  acoustic  pressure  as  a  specifieo  function  of  depth  and 
azimutn  at  some  near-source  radial  distance.  Further,  we  assume  boundary 
conditions  are  specified  at  the  ocean  bottom  ana  surface  ano  for  the  azimuthal 
region  of  interest.  For  brevity  we  do  not  write  these  conditions  iri  the 
following  development,  but  they  are  easily  incorporated  once  physical  models 
for  the  boundaries  are  specifieo.  The  only  spatial  conoitioris  which  we 
explicitly  employ  are  the  obvious  ones  of  boundea  pressure  for  all  ranges  ano 
of  only  an  outgoing  wave  at  large  ranges. 
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With  some  cnoice  of  tne  reference  sounu  speea  cQ<  we  assume  that 
n2(r,  o,  z)  can  be  written  in  the  form 

n2(r,»,z)  =  1  +  en(r*,»*,z*)  ,  (40) 

in  wmch  n  is  a  function  of  oraer-of-magnituae  uriityl2  and  the  scaled 
variables  with  asterisks  are  nondimensional .  The  quantity  n  could  be  taken  as 
the  maximum  relative  deviation  of  c  from  cQ  *nicn  typically  is  no  more  than 
about  10-2  in  ocean  applications.  The  scaled  variables  are  defined  by 

r*  =  e  kQr,  z*  =  e1'2  kQz.  o*  =  0  E~1/2  o  .  (41) 

Tne  first  two  variables  of  Eq.  (41)  are  chosen  following  Tappert1  who  also 
provides  a  justification  for  the  atove  definition  of  e.  (Other  definitions, 
sucn  as  a  cnaririel  aspect  ratio,  may  be  appropriate  for  certain  propagation 
conditions,  as,  for  example,  in  an  isospeeo  channel  where  no  definition  is 
ideal  in  all  circumstances.)  The  third  variable  contains  ati  ordering 
parameter  a  that  serves  to  account  for  the  rapidity  of  the  azimuthal  variation 
in  n2  (r,  e,  z;  ano,  more  generally,  in  the  solution.  Thus,  it  a  is  of 
order  unity  with  respect  to  e  [denoted  by  a  =  0(1)],  then  the  dimensional 
azimutrial  derivative  r-*  ap/s©  is  comparable  to  ap/az,  so  that  one  or  more 
of  the  three  mechanisms  mentioned  previously  are  accounting  for  substantial 
azimuthal  variations.  If  on  the  other  hand  a  =  0(eW2)f  then  r~l  ap/3©  is 
comparable  to  ap/ar,  arid  azimuthal  variations  are  relatively  smaller.  Our 
development  proceeds  with  a,  or  order  one,  and  extends  to  other  cases. 
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Employing  Eqs.  (40)  ana  (41)  in  Eq.  (1)  and  dropping  the  asterisks 
hencefortn  (so  that  (r,e,z)  now  represent  scaled  variables),  for  p  =  p(r,e.z), 
we  obtain 


9  1  * 

e  (Prr  +  -p  pr)  +  e (pZ2  +  P9e)  +  (1  +  gu)  P  =  0  ,  (42) 

where  subscripts  denote  partial  derivatives  with  respect  to  the  scaled 
variables.  Motivated  by  Eqs.  (2)  ana  (6),  we  apply  the  method  of  multiple 
scales  by  seeking  a  solution  of  Eq.  (42)  in  the  form 


Plr ,e,z)  =  P(p,r,o,z;e)  , 


where  p  =  r/e.  With  Eq.  (43),  Eq.  (42)  can  be  written  as 

2 

fp  +  pi  +  e  r  2P  +  -  p  +  P  +^ttP  +nP]  +  e2[P  +  —  P  j  b  0  .  (44) 
lpp  rj  eLt  rp  r  p  zz  ”2  e©  n  J  1  rr  r  rJ  '  ' 


We  note  that  our  results  are,  in  fact,  unchanged  if  (for  example)  the  term 

er-^  is  written  as  p~^  Pp,  but  the  analysis  woula  be  more 
involved.  We  next  assume  an  asymptotic  expansion  of  P,  i.e., 

oo 

p  /  y  6n  P(rl)(p,r,©,z)  ,  e  ->  0  ,  (4b) 

n=0 

ana  inserting  Eq.  (4b)  iri  Eq.  (44)  produces  a  sequence  of  equations.  Tne 
first  three  of  which  are 

^p(O)  ^  p(0)  +  p(0)  =  0  (46) 

PP 
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Equation  (b2),  with  a  =  1  Lwhicn  can  be  chusen  without  loss  of  generality  when 
a  =  0(1) j,  is  precisely  Eq.  (9),  the  fundamental  narrow-angle  3D  PE  of  Ref. 

1.  With  tne  definition 
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p(!J  =  A(i)  eip  =  u(1)(r,G,z)  e1p/r1/2 


(53) 


We  insert  Eos.  (45),  (51),  ana  (53)  into  Eq.  (46)  ana  require  P(2)  to  be 
bounaea  tor  all  p.  It  follows  as  before  that  uU)  must  satisfy 


—2 1 U 


(1)  =  u ( 1 ) 


z  z 


2 

a 

7 


,0) 

0© 


+  nU 


U)  ♦  ut0) 


rr 


4r^ 


(5< 


Our  results  thus  far  are  summarizea  by 

p  ~  ~j2  ^  +  eUUj  +  0(e2)j  ,  £+0  ,  (5 

wnere  u(°)  ana  Ll(^)  are  obtainea  from  Eqs.  (52)  ana  (54),  respectively, 
tor  a  of  oraer  one.  When  a  =  0(e^^),  the  secono  term  on  the  right  of 
Eq.  (52)  Luna  Eq.  (54) j  is  aosent  aria  is  replaced  by  r~2  U^g). 


Having  obtained  an  analogue  of  the  narrow-angle  3D  PE,  Eq.  (9),  we  next 
ootain  a  wider  angle  version  corresponaing  to  Eq.  (20).  Differentiating 
Eq.  (52)  with  respect  to  r  yields 


,(0) 

rr 


i 

7 


,{0) 

zzr 


2a2 

7" 


j(0) 

80 


+  4u(°> 

7  eer 


|(0) 


nU 


(0) 


(56 


It  we  aefine  an  operator 


n  = 


ar 


(57) 


ana  use  Eqs.  (56)  ana  (57)  in  Eq.  (54),  we  fine 
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Adding  Eq.  (52)  and  e  times  Eq.  (58)  gives 

2 

-2UJr  =  Uzz  +  \  U«e  +  nU  +  cnU  "  £^u(i)  >  (&9) 

where  we  have  defined 

U  =  U(0)  +  eU^  .  (60) 

The  last  term  on  the  right  of  Eq.  (55)  should  be  droppeo,  as  is  consistent 
without  neglect  of  0(e2)  terms.  With  this  omission  ana  using  Eq.  (57), 

Eq.  (59)  becomes 


The  resulting  wide-angle  3D  PE  from  our  asymptotic  derivation  is  given  by 
Eq.  (bl)  when  a  is  oraer  one.  Using  Eqs.  (40)  ano  (41),  it  follows  that 
Eq.  (61)  with  a=l  agrees  exactly  with  Eq.  (20)  but  with  the  additional  last 
three  terms  on  its  right  side.  Tnese  tnree  terms  are  multiplied  by  c 
(typically  about  lO-^)  ana  ao  not  involve  a  radial  derivative  of  the 
solution  so  they  are  in  some  sense  less  significant  than  the  remaining  terms 
in  Eq.  1 6 1 ) .  In  fact,  tney  can  be  identified  precisely  with  contributions 
that  were  argueo  as  small  in  the  derivation  of  Eq.  (20).  Specifically,  the 
first  two  of  those  terms  correspond  to  those  neglected  via  Eqs.  (37)  ariu  (38) 
in  the  factorization  approximation.  Furthermore,  the  last  term  corresponds  to 
those  dropped  through  Eqs.  (24)  ana  (25)  in  the  farfieia  approximation.  In 
this  way,  it  follows  that  Eq.  (20)  is  an  apparently  consistent  approximation 
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or  physical  interest  to  the  class  of  3D  PEs  representeo  by  Eq.  (61).  We  note 
tnat  trie  result  of  tne  analysis  for  a  =  0(el/2),  i.e.,  for  relatively  weak 
azimuthal  variations,  may  be  seen  from  Eq.  (61)  by  setting  a2  =  E  ana 
aropping  the  two  0(e2)  terms.  Estimation  of  the  appropriate  magnituae  for  a 
in  any  specific  application  depends  on  detailed  consideration  of  the  three 
azimutna I ly-directive  mechanisms  mentionea  previously. 

A  VALIDITY  TEST 


Trie  accuracy  of  the  3D  wiae  angle  PE  has  been  examinee  by  Schultz,  Lee, 
ana  Jackson5  using  an  exact  solution  test.  Their  exact  solution  u(r,p,z)  is 
required  to  take  the  form 


u(r,o,z)  =  sin(nz)  eime  0(r)  , 

where  p(r)  satisfied  the  differential  equation 


(63) 


dp 

or 


^  ko  m2/(ko r^\ 


(63) 


For  appropriate  choices  of  z  =  an  integer  multiple  of  u,  m=l,  ana  a 
solution  of  Eq.  (63),  the  expression  of  (62)  satisifes  the  3D  wide  angle  PE 
equation  (20).  On  tne  other  hanu,  Eq.  (20)  is  solved  by  an  implicit  finite 
difference  method  that  discretizes  Eq.  (20)  into  a  large  sparse  system  of 
equations.  This  system  was  solved  by  a  Yale  sparse  tecnmque  wnose  solution 
compares  favorably  with  the  exact  solution  (62)  at  every  range  and  at  every  p 
sector.  Results  of  comparison  have  been  reported  in  Ret.  5. 
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We  want  to  establish  the  validity  of  Eq.  (20)  and  compare  the  solution  of 
Eq.  (20)  witn  a  known  reference  solution  for  an  application  (as  reported  by 
Baer?).  Itns  application  problem  deals  witn  a  profile  tnat  can  be 
calculated  by  the  formula  (c(r,©,z)  =  cm(2i  +  (0.001)  rsin(»),  wnere  cm(z) 
takes  on  the  values  described  by  the  table  below  in  the  vertical  plane  at  0°. 


z  (m) 

c(z)  (m/s) 

0.0 

1536.5 

15F.4 

1539.243  ~ 

5CFT3 

nOi.ITU 

1015.9 

1471.882 

5587.91  “ 

1549.606 

5HB775T 

lbbb.S’2fc 

In  the  calculation,  the  source  is  placed  at  254  m  below  the  surface  with 
a  frequency  of  2b  Hz,  ana  the  receiver  is  placea  at  815  m.  The  propagation  is 
carried  out  up  to  140  km  in  range.  Results  are  prooucea  in  azimuthal 
sectors.  The  sector  boundaries  are  assumed  absorbing.  Of  particular  interest 
is  the  result  taken  from  the  sector  [-20°,  20° j  at  range  120  km.  Reference 
results  were  reported  by  Baer^  in  tne  same  sector  using  the  split-step 
Fourier  algorithm.  In  figure  1,  the  solid  line  is  the  3D  PE  result,  the 
dasneu  line  is  the  Nx2D  result,  ana  tne  dotted  line  is  the  3D  wide  angle  PE 
result,  which  was  calculated  by  the  Yale  sparse  technique.  It  is  seen  that 
the  3D  wide  angle  PE  result  compares  closely  witn  reference  solutions. 
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RANGE  =  120  km 


Figure  1.  Results  from  Sector  [-20°,  20s]  at  120  km 


CONCLUSIONS 


A  wide  angle  partial  differential  equation  Has  been  developed  to  predict 
the  underwater  sound  propagation  in  three  dimensions.  This  partial 
differential  equation  is  of  the  third  order  in  theory.  It  is  named  after  the 
313  wide  angle  PE  because  the  small  30  PE  is  a  special  case.  The  entire 
development  was  based  on  an  operator  factorization  whose  theory  was  fully 
justified  by  the  operator  analysis  and  supported  by  the  multiple  scale 
analysis.  The  most  Important  resuit  is  the  information  to  indicate  when  and 
where  the  three-dimensional  problem  can  be  solved  two-dimensional ly.  The 
mathematical  validity  was  established  by  Schultz,  Lee,  and  Jackson^  in  their 
numerical  solution;  however,  the  simulated  example  demonstrated  further  the  30 
wide  angle  PE  capability.  This  3D  Wide  Angle  PE  is,  by  far,  a  more  general 
purpose  model  with  useful  flexible  capabilities. 
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3.  DERIVATION,  CONSISTENCY,  AND  STABILITY  OF  AN 
IMPLICIT  FINITE  DIFFERENCE  SCHEME 


Donald  F.  St.  Mary 
University  of  Massachusetts 

Ding  Lee 

Naval  Underwater  Systems  Center 


ABSTRACT:  Parabolic  equation  (PE)  approximations  to  the  reduced 

wave  equation  (Helmholtz  equation)  are  used  extensively  In  the 
prediction  of  long-range  sound  propagation  In  ocean  environments. 
In  two  dimensions  parabolic  approximating  partial  differential 
equations  have  been  traditionally  solved  numerically  via  a  Green's 
function  approach  (Fast  Field  Program)  and  a  Fast  Fourier  Transform 
(split-step).  Recently,  Lee  et  al.  created  an  Implicit  finite 
difference  ( IFD)  program  to  solve  more  general  two-dimensional  PE 
approximations  (those  that  accommodate  wider  angles  of  propagation). 

In  this  paper,  we  present  a  three-dimensional  PE  (encompassing 
small  and  wide  angles)  that  Is  a  third  order  partial  differential 
equation,  and  derive  an  IFD  scheme  to  solve  It  numerically.  The 
numerical  scheme  Is  presented  In  several  different  ocean 
environments,  a  wedge  shaped  region  with  absorbent  bottom  and 
sides,  the  same  region  with  hard  bottom,  and  a  full  360° 
propagating  region  with  soft/hard  bottom.  Matrix  formulations  are 
carefully  worked  out  In  anticipation  of  the  implementation.  We 
derive  the  consistency  of  the  scheme  with  the  original  partial 
differential  equation  and  show  that  the  scheme  Is  second  order 
accurate.  Finally,  we  present  a  discussion  of  the  stability 
properties  that  might  be  exhibited  by  the  scheme. 
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INTRODUCTION 

In  this  section  we  shall  aiscuss  the  derivation  of  the  implicit  finite 
difference  scheme  associated  with  the  wide  angle  three-dimensional  parabolic 
approximation.  Me  shall  prove  that  the  difference  equation  is  unconditionally 
consistent  with  the  partial  differential  equation  and  investigate  the 
stability  of  the  scheme.  The  finite  difference  approximation  is  a 
Crank-Nicolson  type  scheme.  We  shall  show  that  it  has  consistency  properties 
that  are  very  much  like  those  of  the  classical  Crank-Nicolson  scheme  when 
applied  to  the  canonical  heat  equation.  In  this  regard  we  remark  that  the 
straight  forward  explicit  difference  scheme  is  stable  under  certain  conditions 
on  the  parameters  when  applied  to  the  heat  equation,  but  is  unstable  for  all 
combinations  of  the  parameters  even  when  applied  to  the  simplest  of  our  two 
dimensional  parabolic  partial  differential  equations. 

A  CRANK-NICOLSON  TYPE  APPROXIMATION  SCHEME 

The  wide  angle  approximating  parabolic  equation  (PE)  is  given  by 

(1  +  qx  Lj  +  q2  l2)  ~  -  ikQ((P1  -  q^  Lj  +  (p2  -  q2)  l2)u, 

(1) 

LjU  =  [{n2(r,z,«)  -  1)  +  (l/k2)(a2/3z2)]u  ,  L2u  -  [(l/k|r2  )(a2/oo2 )]  u 

u  =  u(r ,z,e) ,  p1  ±  qp  p2  *  q2-  We  shall  let  u"  -  u(rn,zm,o  ),  where 

+  nk,  r  >  0,  zm  «  mh,  q  •  +  fcd,  thus  Ar  *  k,  az  »  h.  Ad  =  d,  the 

no’o  m  fco  ’ 
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limits  of  the  indexes  are  m  =  0,1,...,  M  or  M  +  1,  l  =  0,1,...,  L  or  L  +  i. 
n  =  0,1,....  In  an  abuse  of  the  notation,  we  shall  use  n  in  two  different 
ways:  as  the  index  of  refraction  and  as  a  counter  of  the  number  of  range 
steps.  (The  context  should  make  it  clear  which  is  intended  in  each  case.)  We 
wish  to  derive  a  Crank-Ni colson  type  approximation  to  (1). 

A  standard  way  in  which  the  Crank-Nicolson  approximation  is  derived  for 
traditional  parabolic  partial  differential  equations  is  to  take  the  average  of 
the  classic  explicit  (forward)  difference  approximation  and  the  implicit 
(backward)  approximation.  In  order  to  motivate  tne  application  of  this 
procedure  to  (1),  we  shall  briefly  describe  its  application  to  a  PE  in 
standard  form,  namely,  the  small  angle  Tappert  equation 


Ur  •  CU  +  du2z  ,  c  =  ikQ(no  -  1  )/2,  d  »  i/2kQ 


(2) 


Consider  the  two  stencils  (for  this  demonstration,  u  »  u(r,z)) 


m-1 


m+1 


fUdl 

| 

B 

|| 

fill 

■mi 

‘m-1 


‘m+1 


i 

V 

C 

n  n+1  n  rn+l 

The  first  of  these  is  used  to  make  the  forward  approximation  based  at  the 
point  (rnJ2m);  the  second  is  used  to  make  the  backward  approximation  based 
at  (rn+l'zm)'  The  difference  equations  are 


,,,n+l  ,,n*,,  n  n  .  n  „  n  n  ...2 

(u„  -  uj/k  =  cm  um  +  d(um,.  -  2um  +  um  t  )/h 

in  nr  mm  '  m+i  m  m-1' 


(3a) 


ana 


,  n+1  ,.n+l  ..n+l  .  wr.,n+l  o  n+l  .  n+1,  ,,2 

^um  -  um)/k  “  cm  um  +  d^m+l  “  2um  +  VlJ/h 


3-4 


(3b) 


TP  '’145 

Note  that  the  left-hand  sides  of  these  equations  are  the  same.  The 
Crank-Nicolson  approximation  to  (2}  is  obtained  in  taking  ((3a)  +  (3b ) ) / 2 . 


In  order  to  begin  to  carry  out  this  development  for  (1)  we  need  to  define 


2  2 
U  i  j  1  J  y 

the  forward  and  backward  (in  r)  discretizations  of  — *■  —  and  — *•  (— ). 

azz  9r 


77  'ar' 


In  each  case  we  shall  take  the  centered  difference  in  the  second  order 
variable  and  the  standard  oifference  in  r.  The  key  to  taking  the  forward  ana 
backward  differences  in  r  is  to  keep  the  base  point  of  the  stencils,  the  point 
at  which  the  approximations  are  being  made,  clearly  in  mind. 


The  two  stencils  associated  with  the  z-derivative,  base  point  encircled 
are 


'Wl-V 


'Wi-V 

(forward) 


^rT 


<rn.l-W 


<W*«> 


(rn+l’zm-l,en) 

<rn+l-W 

(backward) 


The  two  difference  approximations  to  urzz  are  equal,  as  above  (the  forward 
and  backward  approximations  to  the  full  differential  equation  (1)  are  not 
equal  though),  ana  have  the  common  value 


Ki.,  -  2C  *  -  «Ci.,  - 

henceforth  we  shall  use  the  central  difference  operator  notation 


(6mu)ii  =  2um,*  +  Um-I ,i  ^ 

A  completely  analogous  development  takes  place  for  the  second  derivative  in 
e.  The  corresponding  stencils  in  the  r,e  directions  are  as  follows: 
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It  is  not  difficult  to  prove  that  for  arbitrary  sufficiently  differentiable 
functions  0(r,z,e) , 


"  2  0+1  2  "1  /  2 

<*rzz>  *  -  (fimuWkh 

if 

"(*rr“C« {m) '  (4) 

-(*rr4.)  (kh2/24)  *  0(k2  *  n2) 

m,  i 

as  h  *  0,  k  »  0  independently  of  the  manner  in  which  h,k  approach  zero.  A 
completely  analogous  formula  holds  for  *>ree,  which  we  shall  use  in  the 
ensuing  development. 


We  shall  now  obtain  the  desired  difference  scheme  by  producing  the 
analogue  of  (3a)  and  (3b)  for  equation  (1).  For  the  first  of  these,  the  base 
point  is  (rnt2m>0]i)  and,  thus,  the  forward  approximation  is  given  by 
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(1  ♦  qjUn2)"  -  1)  (u£  -  ujtt)/k)  *  «liU/k|)  [(«£<0"*1  ' 


*  q2ti/fk0rn)2  [(4“)"  1  - 

*-  m  mJ 


-  1k0lPi  -  IjMln2)"  -  *  (i/k0)(P,  - 

m  i « 

+(i/K  r^)(P2  -  q2)(6?u)  IQ2 

*■  m 

The  backward  approximation  has  base  point  (rn+i ,Zm,oa)  and  is  given  by 
p  n+1  _+i  9  ,  n+l  0  n 

(1  +  qi((n  }m.,  - 1}  <i  -  u^m  +  qi(1/ko}  -  <4u>*  i ^ 


p  „  n+l  •?  n  ^ 

+  q2U/(k0r^+i)  (4fU)  -  («>)  / kd  ) 

L  x'  m  *■  m 


p  n+l  n+i  o  n+l  p 

-  1K0CPX  -  q^K"  )  -  Du,,  ♦  d/k0)(Pi  ~  /h 

m  ^ 

+(i/korn+l)(p2  -  q2)(4u)  I* 

m 

Finally,  the  average  of  these  two  yields  the  Crank-Nicolson  difference 
approximation  system.  After  considerable  simplification  the  system  can  be 
expressed  as  follows: 


(S'"2 >“£!.,  *  inVKS-i  *  KU  - (2,h2)F 


/ p  , h2  vT-rh ,  n+l  .  ,r-,,2  ,  n+l  ,  #Trrn ,  2.  n+l 
-  <2/d  )bl  K,i  (b/h  )um+l,*  +  <bl  /Q  )uma+l 

=  W^Cl.i  +  (b°n/d2>V<ul  +  -  (2/h2)b 

-(2/d2)bOn)u^  +  (b/h2)ud+1>K+  (bOn/a2)ud?a+1 
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where 


b  =  b(K)  =  q,/k*  +  ik(p1  -  )/2k 


2k^  rL  (r  +  k) 


2k0ir  +  k) 


Qo  i  1  ik(p-  -  qo) 

bO  -  bO(r;k)  =.  — f  (±»  + - ± - 2)  + - - - *  (7 

2k£  (r  +  kr  2kQ  rd 

al  =  al (r  ,z,o;k)  -  (1  +  q^  (n2  (r,z,e)  ♦  n2(r+k,z,e) )/2  -  1].) 

+  ikkQ(p1  -  q1)(n2(r+k,z,e)  -  l)/2 

aO  -  aO(r ,z,e;k)  -  (1  +  qj[ (n2(r,z,o)  +  n2(r+k,z,e) )/2  -  1]) 

+  ikkQ(p1  -  q1)(n2(r,z,e)  -  l)/2 

The  bar  over  an  expression  indicates  the  taking  of  the  complex  conjugate. 
Note  also  that  al  ano  aO  are  equal  if  n(r,z,e)  is  independent  of  the  range 
variable  r. 


BOUNDARY  CONDITIONS  AND  MATRIX  FORMULATION 


We  wish  to  express  system  (6)  in  a  convenient  matrix  formulation,  out  the 
precise  form  of  the  coefficient  matrices  aepenas  on  the  boundary  conditions 
imposed  on  the  original  problem.  Throughout  the  discussion  we  shall  assume 
the  surface  bounaary  conditions  u(r,Q,e)  =  0.  Another  standing  assumption  is 
the  initial  condition,  namely,  for  a  given  function  f(z,e),  u(r Q,z,e)  = 


A  frequently  imposed  bottom  boundary  condition  is  u(r,z^+^te)  »  0  (in 
this  formulation  the  bottom  is  at  z  *  2|*|+}),  the  general  assumption 
associated  with  this  condition  is  that  an  artificially  imposed  absorbing  layer 
below  the  ocean  floor  prevents  energy  from  entering  the  water  column. 

Similarly  we  can  consider  propagation  taking  place  in  a  cylindrical  sector 
(pie-shapeo  region)  between  two  azimuthal  angles  denoted  by  eQ  ana  eL+1 
with  an  absorption  region  on  each  vertical  side  of  the  sector.  Thus  in 
aaaition  to  the  conditions 

u(r,zQ,«)  -  0  ,  u(r,zM+1,©)  .  0  ,  (8a) 

we  have 

u(r »z,Oq)  -  0  ,  u(r,z,oL+1)  .  0  ,  (8b) 

for  all  r  >  0,  z,  Oq  <  9  <  This  is  the  case  considered  by  Baer  and 

Perkins  for  small  angle  PE.  Finally,  one  assumes  a  given  sound  profile  at  a 
distance  from  the  source 

u(r0,z,g)  =  f ( z ,8 )  .  (8c) 

The  system  (6)  in  conjunction  with  the  boundary  conditions  (8)  can  be 
expressed  in  a  particularly  nice  symmetric  block  tridiagonal  matrix  form. 
Namely 
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where  each  block  is  M  x  M,  the  diagonal  blocks  are  tridiagonal  matrices  ana 
the  oft  diagonal  blocics  are  diagonal  matrices  with 
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ru 


*U) 


n 

l.» 
J2  ,1 


JM,*J 


(9d) 


The  system  (9)  shall  be  referred  to  symbolically  as 

ar"  un+1  -  aOnun 

We  shall  now  turn  to  more  general  bounaary  conditions  ana  consider  a 
cylindrical ly  shaped  region  0  <  u  <_  2«.  Again  we  shall  retain  the  pressure 
release  top  surface  boundary  condition 


u(r,zQ,«)  -  0  ,  r  >  0  (10a) 

The  bottom  boundary  conditions  are  artificially  located  far  below  the  actual 
bottom  of  the  wave  guide,  but  in  the  current  case  we  assume  that  the  position 

2  *  is  the  actual  interface  between  ocean  floor  and  water  ana  allow  for 
the  possibility  of  reflection  of  rays.  Fbr  given  real  constants  o,  bq,  y 
the  condition  is  given  by 


«u(r,zM,o)  +  snu7(r,zM,o)  =  Y 


B0  *  0 


(10b) 


The  case  of  general  interest  is  a  -  0,  bq  .  1,  T  .  0.  Finally,  we  impose  a 
continuity  condition  on  the  motion  in  the  e  variable,  namely 


u(r,z,0)  =  u(r,z,2n)  ,  uQ(r,z,0)  =  ue(r,z,2„)  ,  r  >  0 


(10c) 
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We  shall  assume  that  eQ  =*  0  and  o^+1  ■  2*. 

First  consider  the  discretization  of  the  boundary  condition  (10b).  We 
shall  use  a  centered  difference  to  approximate  the  derivative  so  as  to 
maintain  the  second  order  character  of  the  approximations  (as  shall  be  seen  in 
the  ensuing  development,  (6)  is  a  second  order  scheme).  We  use 

uz(r,zM,o)  =  [u(r,zH+1,o)  -  u(r,zM_1,e)]/2h 
and  thus  for  (10b)  we  write 

aUM,ji  +  ®o^uM+l,s,  “  uM-l,J/2h  "  T  ’  *0*°  ‘ 

The  term  uJJ+1  %  is  ficticious  in  this  context  (recall  the  bottom  is  at  u”  ^ ) 
ano  can  be  expressed  in  terms  of  the  real  unknowns  ujjj  z  SL  using  (11). 

In  order  to  encompass  (11)  into  the  matrix  formulation  of  the  problem,  set  m  > 

M  in  (6)  ano  make  the  substitution,  from  (11), 

uS+l,z“uS-l,z-2ahuS,^o  +  2h^Bo  • 

to  obtain 

2<“'h2»us'-i.t  *  cl-i  *  <*C  -  (2/hZ> 5  + «w8o> 

-  <2'“2>  *  c  *  ST1'”2 

.  2(b/h2)uJ_1>l  *  (bOn/d2)uJJ^_j  *  (aOj_t  -  (2/h2)  1.(1  *  .h/»0)) 

-  (2/d2)  bo")uJjt  *  (boV)^,.,  *  gh>k  ,  (12) 

9h,k  =  (2y/e0h) (b  -~o  =  2vki(Pi  -  q1)/B0hkQ) 
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Condition  (12),  compared  with  (6),  forces  a  change  in  the  last  row  of  each  of 
the  diagonal  blocks  aln,  aOn,  a-  1,2,..,,L,  of  the  system  (9). 

In  particular  if  one  should  use  a  bottom  reflecting  condition  (11)  in  a 
sector  with  absorbing  vertical  sides  these  row  changes  plus  the  addition  of 

the  g^  ^  vector  would  be  the  only  change  in  (9).  The  resulting  matrix 
system,  which  we  choose  to  denote  by 


ar'n  un+1  -  aO  n  un  ♦  g 


(13a) 


has  the  same  coefficients  as  (9)  except  that  for  the  AOn  we  substitute 


AO 


•n 


aO 


I,* 


where 


aO 


2, A 


20 


•'M-l  ,z 


aO 


•n 

M,t  -1 


1,2 . L,  (13b) 


-0^  »  aOjjJt  -  (2/h2)  b(l  +  *h/eo  -  (2/d2)  bOn 
The  analogous  change  in  constructing  the  Al  from  (12)  is  made, 


I n 

lM,a 


(2/h2)  b(l  +  ah/so)  -  (2/d2)  bln 


and  g  is  the  M  x  L  -  column  vector 

9  =  0 » •  *  •  ,0  *  9^  ^  ]  t  •••»  LO, .  -  -  ,°»9h,k^] 
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The  system  (13),  which  does  not  exhibit  the  symmetric  character  of  (9),  will 
be  shown  to  be  equivalent  to  a  symmetric  system. 

Returning  to  the  main  consideration  of  this  subsection,  namely,  a 
cylinarically  shaped  region  with  bottom  reflecting  condition  (11),  the 
conditions  in  (10c)  can  be  represented  as  follows: 


um,0  *  um,L+l  »  tum,l  "  um,0^d  “  ^um,L+l  "  um,L^d  * 

The  dummy  index  0  ■  l+l  is  only  used  to  help  indicate  direction  of  approach  to 

the  vertical  plane  o  -  0,  i.e.,  u”  «  uj^  Q.  The  two  relations  in  (14) 
reduce  to 


m,0 


<1  +  <,l  m 


(15) 


One  now  re-examines  (6)  in  the  critical  cases  t  -  1,  a  ■  L; 
m  -  1,2,. ..,M.  The  new  system  of  equations  differs  from  (13a)  only  in  the 
first  ana  last  rows  of  blocxs.  It  is  an  (ML)  dimensional  square  system 
block-tridiagonal-like  in  form,  except  for  the  first  and  last  rows  of  blocks 
each  of  which  has  one  additional  block;  i.e., 


where 


ATn  un+1 


AO*'11  un  +  g 


(16) 
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AKi"  BKn  ( 1 / 2 ) BKn  “ 

6Kn  AKjH  Bn 

•  •  • 

W'n  ~  .  .  K  *  0,1 

l 

•  •  • 

n  11  n  n 

BK  AKl_x  BKn 

(1/2) BKn  BKn  AK^n 

BOn,  Bln,  6  are  as  in  (9), 


where,  for 

i  *  1,L  ;  aK^-  aKj  ^  -  (2/h2)  b  -  (3/2d2)  bkn  ,  m  =  1,2 . 

aicj^-  aKM"j_  -  (2/h2)  b(l  +  ah/B0)  -  (3/2d2 )  bKn 

and  for 

-  ^ . L-l  ;  a^,.  -  <>t  -  (2/h2)  b  -  (2/d2)  bl<"  ,  . 

-  W''2)  "d  *  *h'» 0>  -  d/"2)  bk 
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K  -  0,1.  Of  course  since  everything  on  the  right-hand  side  of  (16)  is  known, 
one  could  write  it  (Out  not  the  left-hand  side)  without  the  addition  of  the 
corner  blocks,  using  (15)  directly  ano  making  the  corresponding  changes  in  the 
vector  g. 


CONSISTENCY  OF  IMPLICIT  FINITE  DIFFERENCE  ( IFD) 

A  difference  equation  approximation  to  a  partial  differential  equation  is 
said  to  be  (unconditionally)  consistent  with  the  differential  equation  if  the 
difference  equation  approaches  the  differential  equation  as  the  mesh  size 
approaches  zero,  independently  of  the  manner  in  which  the  mesh  size  approaches 
zero.  More  precisely,'  (6)  is  consistent  with  (1)  if 


n 

b  An+1  +  bT  j.n+1  .  /— r  2  t-  2  rr  >  ,n+ 

77  Vl.t  -J-  V-l  "  I?  b  "  ~J  W  >  V 


n+1 

l 


n  — 

+  j.n+1  *  "5T  An+1  b  .n  bO  .n 

^7  Vl.fi  ~^T  Va+i  ~  ^  Vl,a  "  *hm-1 


K.i  <17> 


bOn  ,n 

2"  ~-J  ™  1  ~  ^  *Vl,t  '  7"  ^m,2.+l 


2  h  2  hnn,  *n  b  ,n 
-  b  -  “7  b0  )  o  -  —  K 
h‘  h 


(1  +  qx  Lj  +  q2  L2)  t>r  -  ikQ  ((pj  -  q1)L1  +  (pg  -  q2)L2  (6 

approaches  zero  as  h,  k  >  0,  independently  of  the  manner  in  which  h,  k 
approach  zero  for  arbitrary  net  functions  p(r,z,e)  having  sufficient 
differentiability.  The  factor  1/k  is  present  since  in  the  derivation  of  (6) 
we  previously  cleared  the  k  from  the  denominator.  In  order  to  help  simplify 
(17)  we  shall  express,  aO,  al,  bO,  bl,  and  b  in  terms  of  their  consistent 
parts.  Let 
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t  . 


Ra  *  Ra(r,z,«;k)  -  1  +  q1  [(n2(r,a,e)  +  n2(r  +  k,z,e))/2  -  l] 

la  -  Ia(r,z,«)  *  ikQ  (p1  -  ^ ) (n2 (r,z,o)  -  l)/2 

RbN  -  RbN(r;k )  .  q2((l/r2)  +  (1 / (r  +  k)2))/2k2 

IbN  -  IbN(r)  -  1(p2  -  q2)/2kQr2 

Rb  -  qx/k2  ,  ID  -  i  (Px  -  qj  / 2kQ 

c  -  c(r ,z,e;k)  .  q1[n2(r  +  k,z,e)  -  n2(r,z,e)]/2  , 

then 

aO(r,z,o;K)  «  Ra(r,z,o;k)  +  kla(r,z,®)  , 
al(r,z,o;k)  -  Ra(r,z,e;k)  +  kla(r+k,z,e)  , 
bO(r;k)  -  RbN(r;k)  +  klbN(r)  , 
bl(r;k)  -  RbN(r;k)  +  kIbN(r  +  k) 
and 

b  -  RB  +  klb 
It  follows  that 

(1  +  qx  Lj )  Pr  -  (Ra  -  c)  tr  *  Rb*rzz 
1kQ(Pi  -  qj)  -  2(1  ap  +  Ib*zz) 
ikQ  (P2  -  q2)  L20  «  2IbN*M 

and  that  the  standard  Taylor  approximation  in  the  r  variable  yields 
c(r,z,e;k)  *  q^Jlr.z.e)  k/2  +  0(k2) 
ana 

RbN(r;k)  «  q2/k2r2  -  kq2/k2r3  +  0(k2) 

Now  (17)  can  be  expressed  in  the  form 
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n 

) 

m,j 


rr  (^oo)  +  2 lb(0  )  +  2IbNn  (0oa) 


Each  term  appearing  in  brackets,  [•♦♦]»  can  be  expanded  using  a  standard 
Taylor  approximation,  the  centered  difference  approximation,  or  (4),  thus  we 


obtain 


Ra"  i  \l*rr"  k/2  +  °<k2)l+  (Rb  "  kIb)  krzz>  +  (0rrzz)  k/2 
ni*J  L  rr  m,j  J  L  m.j  m.j 

♦  0(h2  +  k2)l  +  (RbNn  -  klONn)  [(0  )  +  ($rree)  .  */ 2  +  0(d2  +  k2) 

j  L  (r  |  j  j 

♦  k  (ry)n  k/2  +  0(k2)l  (0  )"  -la"  [(0r)  .  k  +  0(k2)l 

L  r  m,j  J  m,j  L  m,j  J 

-  ^[^zzC.J  *  *  [<j  *  °(k)]  k  +  0(k2>] 


-  [(IbN  )"  k  +  0(k2)  (0  )  +  0(k)  +  0(d2)l  -  2IbNn[(0  )  +  0(d2) 

L  r  m.j  m,j  J  L  m.j  J 
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! 


n 

I 


-“bWr22t  .  -  <2 -f?  “rw1  .  *2II><*2Z>  ."2I»N  «„> 

rn,j  kQrn  ,wm,j  m,j 

(k/2)  [Ra^rr  *  mmi  *  RbNJrrw  *  qi(n2)r  *r 

-  (2q2/k2r3)  *r98  -  2  (la4r  *  lb  tm) 

-  2IbN  d  -  2Ia.d  -  2IbN«  )  1*  0(b2  *  k2  *  d2) 

m,  jJ 

(k/2)  3[{1  +  qj  L1  +  q2  L2)*r  -  ikQ  ( (px  -  qL)  x  +  (p2  -  q2 ) L2 )  #]/ar 

J  m,  j 


+  0(h2  +  k2  +  d2)  ,  (18) 

where  the  last  equality  uses  the  fact  that  kctfrr.»  0(k2)-  It  follows 
immediately  from  the  equality  of  (17)  and  (18)  that  in  the  range  dependent 
index  of  refraction  case,  n  *  n(r,z,e),  the  Crar.k-Nicolson  difference  scheme 
(6)  is  unconditionally  consistent  with  the  partial  differential  equation  (1). 


Further,  the  truncation  error  or  local  discretization  error  can  be 
obtained  as  the  magnitude  of  the  difference,  at  a  point  (rn,zm,e  )  between 
the  differential  equation  and  the  difference  equation  both  evaluated  with  the 
net  function  4  =  u,  the  exact  solution  of  the  partial  differential  equation. 
Again,  the  equality  of  (17)  and  (18)  yields  Immediately  the  result  that  the 
local  discretization  error  of  (6)  is  Q(h2  *  k2  +  d2). 

TRANSFORMATION  AND  NONSINGULARITY  OF  THE  DIFFERENCE  SYSTEM 

The  system  (16)  is  a  square  system  of  equations  having  M  x  L  unknowns. 

We  shall  show  that  a  solution  always  exists,  i.e.,  that  al  M  is  nonsingular, 
under  practical  assumptions  on  the  parameters  in  (1).  Looking  at  (16a)  one 


.’■V 


?■  \sl 

"V 

1 

'V- 
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can  see  that  the  last  row  of  each  AKjn  causes  nonsymmetry  to  enter  (16). 
We  shall  show  that  a  very  simple  transformation  of  the  original  system 
replaces  it  witn  an  equivalent  system  in  which  the  AK']n  are  symmetric,  K 

J 

=.  0,1. 


First,  we  shall  rewrite  (16)  in  a  manner  which  brings  out  the  structure 
of  tne  individual  blows  AK"n.  Let  DkJ  be  trie  M  x  M  diagonal 
matrices  having  diagonal  elements  as  follows: 


j  -  diag  [aKj)jt  aK^,  ....  aKj^]  ,  j  -  1,...,L  ;  K  =,  0,1, 


DK 

and  T  be  the  M  x  M  tridiagonal  matrix 

'2  -1 

-1  2  -1 


-1 


2 

-2 


-1 

2(1  ♦  ~) 
* o 


then  for 


j  *  1,  L.  AK.n  -  DK1]  -  (b/h2)T  -  (bKn/d2)  (3/2)1 

J  J 


and  for 


j  «  2,  ....  L  -  1,  AK^n  =  DK"  -  (b/h2 )T  -  (bKn/d2 ) (21 )  ,  K  =  0,1, 

where  I  is  the  M  x  M  identity  matrix.  Let  DKn  ana  J  denote  the  ML  x  ML 
(block)  diagonal  matrices 


DKn  =  diag  [DkJ,  DkJ,  ...,  DKJ  ,  K  «  0,1 
J  =  diag  [I,  l,  ...,  I] 


\ 


) 


(J  is  the  ML  x  ML  identity  matrix],  and  let  S  be  the  block  triaiagonal  matrix 
with  two  additional  blocks 


then 

AK"n  .  DKn  -  (b/h2 )TJ  -  (oKn/d2)I  S  K  *  0,1 

where  we  are  using  some  obvious  scalar  block  multiplication  of  M  x  M  matrices 
ana  ML  x  ML  matrices. 

We  shall  now  pursue  the  transformations  alluded  to  above.  The 
nonsymmetry  obviously  arises  from  the  involvement  of  the  matrix  T.  Let  P  be 
the  M  x  M  diagonal  matrix,  i.e., 

P  -  diag  [1,  ....  1,  1/V?]  , 

then  T  =  P'^SP,  where  S  is  the  M  x  M  symmetric  tridiagonal  matrix  having 
exactly  the  same  entries  as  T  except  in  the  lower  right  2x2  block  where  S  is 
of  the  form 

-sfl 

2(1  +  ah/eQ) 


2 

S(lower  block)  = 

-V? 


Let  P,  SQ  be  the  ML  x  ML  block  diagonal  matrices 

P  =  diag  [P,  ....  PJ  ,  SQ  =  diag  [S,  S] 

then 

AK"n  =  DKn  -  (b/h2)  -1  SPD  -  (bKn/d2) IS 

-  P~]  DKn  -  (b/h2)  P_1S0P  -  (bKn/d2)  P~XSP 

*  P~l[0Kn  -  (b/h2)  SQ  -  (bKn/d2)S]P 

=  P_1  AtK"n  P  ,  K  -  0,1 
Thus  the  system  (16)  may  be  expressed  in  the  symmetric  form 

A^r"n  vn+1  =  Ato"n  vR  +  Pg  ,  (19) 

where 

„n  .  p-V 

and 

AtK"n  =  DKn  -  (b/h2 )SJ  -  (bKn/d2 ) IS  ,  K  a  0,1  .  (19a) 

In  order  to  obtain  the  nonsingularity  of  the  system  we  need  to  derive 

II 

conditions  under  which  the  matrix  in  (19a)  is  nonsingular.  Let  A  a  Aj.1  , 
n  be  fixed  but  arbitrary,  and  suppose  there  exists  an  ML-vector  y> 
y  =  (Yp  Y2»  ti_)»  •••»  M  vectors,  such  that  r  *  0.  then 

y*  y  =»  i  -e. , 


(b/h2)  y*  gY  +  (bln/d2)  y*Sy  =  T*^lnv 

The  quadratic  terms  y*SQy,Sy*  y  are  real  since  SQ,  S  are  real  symmetric 
and  it  is  easy  to  see  that  y*$y  1  0.  Now  separate  the  equation  into  two 
equations  by  taking  its  real  part  and  its  imaginary  part,  then  eliminate 
y*Sqy  from  the  two  equations.  The  result  is  the  expression 
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1/1.1  ,  q2^Pl  "  ql}  /  1  ,  ql(p2  “  q?) 

?  (7  7“}  “77 - 7”  ~77 

rn  rn+l  kod  rn+l  V 


Y*Sy 


L  M 

?  n  2  n+1 

(P,  -  q,)  >  >  (1  +  q,  ((7)  -  (r/)  )/ 2 

4  *—*  m,j  m,  j 

j  «  1  m  -  1 


(m) 


,  (20) 


rj-  Wf . j«)  .  We  wish  to  state  conditions  under  which  it  is 

impossible  for  (20)  to  hold  except  for  y  the  zero  vector.  The  first  condition 
we  shall  state  is  where  the  index  of  refraction  is  slowly  varying  in  range. 
This  is  a  standard  assumption  which  is  frequently  utilized  long  before  this 
point  in  a  development  such  as  this  in  the  general  area  of  underwater 
acoustics.  We  implement  the  condition  here  to  imply  that  the  difference 
involving  n2  in  (20)  is  small.  The  standard  choices  of  the  p,q  parameters 
are  P1  =  p2  *  3/4,  qj^  -  q2  -  1/4;  pj_  -  p2  -  1/2,  qj_  -  q2  -  0, 
or  values  close  to  these.  Under  such  circumstances,  the  right  side  of  (20) 
can  be  seen  to  be  close  to  the  magnitude  of  the  original  vector  y  ,  which  can 
be  taken  to  be  unity  (if  y  ^  0),  and  the  left  side  of  (20)  can  be  made 
arbitrarily  close  to  zero  by  choosing  appropriate  range  step  sizes  k.  Thus, 
we  conclude  that  under  appropriate  conditions  on  the  parameters  that  y  must  be 
the  zero  vector,  i.e.,  the  difference  system  is  nonsingular. 


STABILITY 

We  now  turn  to  the  question  of  the  stability  of  the  scheme.  A  difference 
system  is  said  to  be  stable  if  an  error  (initial,  round-off,  etc.)  made  at  the  nth 
step  does  not,  magnify  uncontrolled  in  its  propagation  to  the  (n  +  l)th.  In  the 
simplest  cases  this  translates  into  a  need  to  show  that  for  a  system  of  the  form 
ul1+1  =  Bun  +  g,  the  eigenvalues  of  B  are  less  than  or  equal  to  unity  is  magnitude. 
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In  the  current  development  the  complicated  nature  of  the  matrix  Bn, 
from  (19),  is 

up  -1 

Bn  -  (7^1  )  (At0"n) 

which  makes  it  very  difficult  to  attack  the  question  of  the  magnitude  of  the 
eigenvalues  of  Bn.  Thus  we  shall  pursue  a  more  heuristic  discussion  of  the 
stability  question.  To  simplify  matters  somewhat  more  we  shall  consider  only 
the  small  angle  PE,  i.e.,  px  «  p2  .  1/2,  q:  -  q2  -  0,  and  we  shall 
assume  the  index  of  refraction  is  constant,  i.e.,  n(r,o,z)  -  n.  Thus  the 
partial  differential  equation  becomes 


ik. 


(n2  -  1)  u  + 

0 


h  92U 

k^r7  a©7 


(21) 


The  von  Neumann  or  Fourier  series  method  of  analyzing  stability  is  a 
method  that  actually  applies  only  to  linear  difference  equations  with  constant 
coefficients  and  then  only  to  initial  value  problems  with  periodic  initial 
data.  In  practice  the  method  is  widely  used  outside  of  this  narrow  band  of 
problems  and  it  frequently  gives  useful  results. 


We  begin  the  method  upon  assuming  that  a  solution  of  the  difference 
equation  (6)  is  given  by 


u 


n 

m. 


ei'w(mh)  e  i  v  ( j d ) 


(22) 


where  5  =  eaK,  o  a  complex  constant.  We  seek  conditions  under  which  (22) 
satisfies  (6)  and  | C  |  <_  1  for  all  n.  Frequently  c  is  called  the 
amplification.  Substituting  (22)  into  (6)  and  simplifying  one  obtains 
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-b  (1  -  cos  u>h )  +  —  (1  -  cos  vh)  -  k2  (n2  -  1)]  +  1 

h  d  rn+l 


-b  (L?  (1  -  cos  wh)  +  7  (1  -  cos  vh)  -  k2  (n2  -  1)]  +  1 

h  d  rn 


(23) 


2  2 

where  we  have  used  the  fact  that  bl  «  b/(r  +  k)  ,  bO  -  b/r  ,  and  al  -  aO  *  1 
+  bk£(n  -  1)  due  to  the  assumptions  enumerated  above. 


Case  1 : 


If  the  inciex  of  refraction  n  -  1  then 


2 


1  +  (ft-) 


— 7j -  (1  —  COS  U)h)  *■  — 5 — tv  (1  -  cos  vh) 

a 


n 


(24) 


1  +  ^ 
o 


^-k  (1  -  cos  wh)  +  "  «2x —  (1  -  cos  vhl 
•h  d  1 


We  do  not  mean  to  imply  that  there  is  no  dependence  on  the  index  n  in  the  left 
sioe  of  (24).  We  are  here  concerned  only  with  the  dependence  of  E  on  n  which 
might  be  of  magnitude  greater  than  unity.  Indeed  separation  of  variables  in 
(21)  indicates  that  its  solution  as  a  function  of  r  has  a  factor  of  the  form 


R(r)  -  e,k°P‘ln  -  1!r  ,*/r  , 

u,x  constants.  It  is  apparent  from  (24)  that  jc|  >  1  which  suggests 
instability  of  the  scheme  in  this  simplest  of  cases.  Extensive  numerical 
computations  have  not  yet  been  carried  out  in  connection  with  the  scheme.  If 
it  should  turn  out  that  the  computations  indicate  no  instability  for  this 
particular  case  then  it  will  be  necessary  to  abandon  the  Fourier  method  and 
analyze  stability  via  a  different  approach. 
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Case  2: 


If  n  >  1  and  k0(n£  -  1)  is  such  that 


X  *  (1  -  cos  uh)  +  \  "y  (1  -  cos  vh) 

d2r2n 


k*  (n2  -  1)]  <  0  , 


*n+l  <  xn  <  and  thus 


1  +  ^  Xn 
1  *  lib*  Xn+1 


Hence  the  method  indicates  stability. 
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4.  APPLICATION  OF  THE  YALE  SPARSE  TECHNIQUE  TO 
SOLVE  THE  THREE-0 IHENSIONAL  PARABOLIC  EQUATION 


Martin  H.  Schultz 
Yale  University 

Ding  Lee 

Naval  Underwater  Systems  Center 

Kenneth  R.  Jackson 
University  of  Toronto 


ABSTRACT:  The  Yale  University  sparse  matrix  technique  is  an  effi¬ 
cient  method  for  solving  large  sparse  systems  of  linear  equations 
such  as  those  that  arise  at  each  step  In  the  numerical  Integration 
of  the  stiff  system  of  ordinary  differential  equations  resulting 
from  the  application  of  the  finite  difference  discretization  to  the 
three-dimensional  parabolic  wave  equation.  We  discuss  the 
procedure  of  a  special  technique,  the  Conjugate  Gradient  method  for 
Normal  Equations  (CGNE)  together  with  Its  advantage  for  solving 
three-dimensional  underwater  acoustic  wave  propagations. 
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INTRODUCTION 


Interest  in  applying  the  parabolic  equation  (PE)  approximation  to  solve 
three-dimensional  ocean  acoustic  wave  propagation  is  on  the  rise.  A 
three-dimensional  parabolic  equation  (3D  PE),  originally  introduced  by 
Tappert,*  dealing  with  small  angle  propagation  was  solved  by  Baer-Perkins2 
effectively.  Baer-Perkins  solved  the  3D  PE  by  means  of  the  split-step 
algorithm  extended  to  three-dimensional  calculations.  Their  efficiency  in 
calculation  is  to  specialize  the  problem  into  N  by  two-dimensional  problems  (N 
x  2D  algorithm).  A  second  three-dimensional  wave  equation  was  recently 
developed  by  Siegmann,  Lee,  and  Kriegsmann3  that  offers  the  wide  angle 
capability  (3D  wide  angle  PE).  Reference  3  showed  that  the  3D  PE  is  a  special 
case  of  the  3D  wide  angle  PE.  Since  the  3D  PE  is  a  special  case  of  the  3D 
wide  angle  PE,  we  proceed  only  to  seek  the  solution  of  the  3D  wide  angle  PE. 
Bayliss,  Goldstein,  and  Turkel4  used  the  sparse  matrix  technique  and 
effective  preconditioning  to  solve  the  Helmholtz  equation.  For  the  solution 
of  our  problem,  we  introduce  the  Yale  University  sparse  matrix  technique.  A 
brief  discussion  on  the  Yale  sparse  technique  will  be  given  in  the  next 
section.  In  order  to  set  up  the  3D  wide  angle  PE  in  the  form  solvable  by  the 
Yale  sparse  technique,  we  apply  an  implicit  finite  difference  scheme  to 
formulate  the  3D  wide  angle  PE  into  a  finite  difference  equation.  Numerical 
solution  to  this  implicit  finite  difference  equation  is  carried  out  by  the 
convergent  Crank-Nicolson  scheme.  A  section  is  devoted  to  discuss  the  finite 
difference  formulation.  To  support  the  validity  of  the  solution,  two  examples 
are  included:  one  demonstrates  the  exact  solution  test  and  the  other  exhibits 
an  application  that  had  been  considered  by  others. 
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AN  APPLICABLE  YALE  SPARSE  MATRIX  TECHNIQUE 


A  linear  system  of  the  form 


Ax  -  f  (1) 

can  be  solved  by  two  classes  of  methods,  the  direct  method  and  the  iterative 
(indirect)  method,  where  A  is  a  square,  nonsingular  matrix  of  order  N,  and  x 
and  f  are  vectors.  All  direct  methods  employ  the  Gaussian  elimination 
procedure,  which  is  very  suitable  for  dense  systems  but  has  limited  usefulness 
for  solving  sparse  systems  because  excessive  memory  storage  is  required  for 
large  N.  This  is  where  the  sparse  matrix  technique  plays  an  important  and 
useful  role  in  obtaining  an  efficient  solution.  These  large  sparse  matrices 
usually  come  from  the  Method  of  Lines  (MOL)  discretization  of  partial 
differential  equations.  There  are  many  techniques  introduced  to  solve  the 
sparse  system  and  an  overview  of  recent  developments  of  these  methods  can  be 
found  In  Ref.  5  (Elman).  Among  these  methods,  a  particular  effective  method 
applicable  to  solve  our  three-dimensional  wide  angle  underwater  acoustic  wave 
equation  is  the  Yale  University  sparse  matrix  technique  package. ^  One  of 
the  sparse  techniques  contained  in  the  package  is  known  as  the  Conjugate 
Gradient  (CG)  method,^’®  which  has  been  developed  to  solve  symmetric, 
positive-definite  systems  iteratively  with  great  efficiency.  In  theory,  these 
iterative  methods  must  converge  and  must  converge  fast  for  efficiency.  In 
practice,  conventional  iterative  methods  require  the  estimate  of  some  kind  of 
parameter  (e.g.,  the  extreme  eigenvalues  of  the  matrix  operator  A)  for  fast 
convergence.  Without  this  estimate  one  has  no  idea  how  fast  his  applicable 
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iterative  technique  converges.  This  is  a  drawback  of  most  Iterative  methods. 

The  CG  method  minimizes  a  certain  norm  in  each  step  and  is  in  a  sense  optimal 
over  a  class  of  iterative  methods.  Since  the  system  is  sparse,  the  operations 
are  inexpensive  and  easy  to  implement.  All  these  properties  make  the  CG  a 
strong  candidate  as  one  of  the  most  robust,  rapid  convergent  iterative 
methods.  This  is  the  reason  we  introduce  it  tp  solve  the  wide  angle 
three-dimensional  partial  differential  equation.  In  application,  the  CG 
method  is  effective  for  solving  symmetric,  postive-definite  problems.  In 
fact,  the  partial  differential  equation  governing  the  ocean  wave  propagation 
with  wide  angle  capability  does  not  always  result  in  a  positive  system.  On 
the  contrary,  it  results  in  a  complex  system.  The  CG  method  cannot  be  used 
unless  an  effective  preconditioning  technique  is  applied.  The  efficiency  of 
the  application  of  the  CG  method  to  solve  the  3D  wide  angle  PE  can  be  enhanced  by 
preconditioning.  These  preconditioning  techniques  solve  the  system 

Ax  =•  f 


by  an  equivalent  system 


Q  2Ax  *  Q_1f  ,  (?) 

where  Q“^  is  in  a  sense  an  approximation  of  A~*  so  that  Eq.  (2)  can  be 
solved  very  economically  because  the  actual  operation  of  Q^A  need  not  be 

performed  explicitly,  and  at  the  same  time  the  condition  number  of  A  is 

improved.  Since  our  resulting  1401  discretization  of  the  3D  wide  angle  PE  is 

neither  a  real  system,  nor  has  the  positive-definiteness  property,  we  use  the 
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A*  for  q_1  as  the  preconditioning  matrix.  We  then  extend  the  technique  to 
handle  a  complex,  nonsymmetr ic  system  whose  solution  is  to  be  shown 
effective.  The  method  we  consider  here  is  recognized  as  the  application  of 
the  CG  method  to  the  normal  equation. 

We  begin  by  dealing  with  the  solution  of  the  system  of  equations  of  the 
form  of  Eq.  (1) ,  i  .e. , 


Ax  *  f  , 

where  A  is  a  nonsingular,  square  matrix  with  complex  elements.  This  problem 
is  equivalent  to  the  normal  equation, 

A*Ax  -  A*f  ,  (3) 

where  A*  is  the  complex  conjugate  of  A.  This  suggests  that  one  natural  way  to 
solve  a  nonsyminetric  system  is  by  applying  preconditioning  to  the  original 
system  and  solving  the  equivalent  system  (3),  provided  no  extra  work  is 
introduced. 

In  theory,  when  the  CG  method  is  applied  to  solve  system  (3),  the  iterate 
x j  minimizes  the  residual  norm.5’  One  member  of  the  CG  family  that  can  be 
used  to  solve  system  (3)  is  known  as  the  Craig's  method'7  and  was  proposed  by 
Hestenes.^  In  this  implementation,  the  iterate  x.  minimizes  the  residual 
norm.  This  is  the  method  we  used  for  our  underwater  applications  and  we 
further  extend  this  application  to  complex  arithmetic. 
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CRAIG'S  ALGORITHM 

The  computation  of  the  Craig's  method  involves  5  steps,  i.e., 

a,j  *  (r.j  ,r^  )/(p^»P-j)  » 

*1+1  -  xi  +  aipi  ’ 
ri+l  “  ri  “  aiApi  ’ 
bi  *  (r>l’r1)/(ri,ri)  • 

Pi+1  -  A*Vl  *  bipi  ’ 

where  rQ  ,  f  _  Ax0,  x0  is  chosen  arbitrarily,  and  p0  -  A*rQ.  The 
above  loop  is  repeated  starting  with  i  =■  0  until  convergence. 

The  work  per  loop  requires  5N  multiplications,  plus  2  matrix-vector  products. 
Besides,  only  4N  storages  are  required  for  the  vectors  x,  r,  p,  and  Ap. 

When  oealing  with  the  solution  of  system  (1),  we  apply  the  preconditioning 
technique  to  transform  system  (1)  into  system  (3).  Then,  Craig's  method  is 
used  to  solve  system  (3).  It  is  natural  to  think  about  the  need  for  explicit 
computation  of  A*A.  The  advantage  of  using  Craig's  method  is  that  the  A*A 
need  not  be  carried  out  explicitly.  This  has  been  clearly  demonstrated  in  the 

computation  procedure. 


THE  3-DIMENSIONAL  WIDE  ANGLE  WAVE  EQUATION 

AS  we  mentioned  In  the  previous  section,  there  exist  two  different  types 
of  three-dimensional  wa.e  epuattons  as  a  result  of  the  PE  approximation.  One 
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is  the  three-dimensional  parabolic  wave  equation  (the  30  PE),  originally  used 
by  Tappert*  to  derive  the  standard  two-dimensional  PE.  Solution  to  the  3D 
PE  has  been  developed  by  Baer  and  Perkins^  using  the  Split-step  algoritnm. 

The  second  type  is  the  three-dimensional  wide  angle  partial  differential 
equation,  developed  by  Siegmann,  Lee,  and  Kriegsmarin.^  (We  refer  to  this 
equation  as  the  30  wide  angle  PE.)  We  chose  to  concentrate  on  the  solution  to 
the  30  wide  angle  PE  because  the  30  PE  is  a  special  case.  We  want  to  remark 
why  we  are  motivated  to  solve  the  30  wide  angle  PE  instead  of  30  PE;  in 


particular  the  application  of  the  Yale  sparse  technique.  In  this  event,  the 
vertical  angle  of  propagation  is  roughly  larger  than  15*,  due  to  the  irregular 
nonzero  boundary  conditions,  or  other  environmental  properties  where  the  fast 
Fourier  transform  (FFT)  is  not  easily  applicable,  this  is  why  a  general 
purpose  solution  is  needed. 


Now,  consider  the  3D  wide  angle  PE. 


—  u  =  -ik  +  ik 
3r  l  o 


i  ♦  Pl  [,W)  - 1  ♦  j  j] .  P; ^ 

°  1  *  dj  n2(r,»,z)  -  1  *  -j  tq  — 

1  k^  ^  1  (k  r)d  3/ 

n  J  n 


ko  92 


(kQr ) 


where  n(r,o,z)  is  the  index  of  refraction  and  kQ  is  the  reference  wavenumber. 


Note  that  when  p1  =.  p2  =,  1/2  and  .  q2  -  0,  Eq.  (4)  reduces 
exactly  to  the  3D  PE.^  Using  the  split-step  algorithm  to  solve  Eq.  (4)  is 
not  easily  applicable.  One  can  easily  see  that  an  alternate  general  purpose 


technique  is  needed  to  solve  Eq.  (4).  One  approach  that  was  considered  was  to 
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multiply  both  sides  of  Eq.  (4)  by  the  operator  in  the  denominator  of  the  riqht- 
hand  side  of  Eq.  (4).  The  following  was  obtained 


1  +  <li 


fn2(r,e,z)  -  1  +4-  -C] 

L  ko  92  J 


1ko  prV 


n  (r,e,z)  -  1  + 


1  3* 

?“? 


(kQr )  as 


3r 


1  321  .  .  _  v  1 

<S^2'(T 


(k0>")  as' 


(5) 


Eq.  (5)  is  not  a  PE,  but  a  third  order  partial  differential  equation  known  as 
the  pseudo-differential  equation.  (A  reminder  to  the  reader  here  is  that  Eq. 
(4)  is  called  the  3D  wide  angle  PE  because  the  3D  PE  is  a  special  case  and  the 
terminology  PE  is  a  very  familiar  term.) 


In  solving  Eq.  (5),  St.  Mary  and  Leey  attempted  to  seek  a  finite 
difference  solution.  Their  analyses  indicate  a  restrictive  stability 
condition.  For  this  reason,  we  attempted  a  similar  implicit  finite  difference 
scheme  as  used  for  the  30  wide  angle  PE  because  of  its  favorable  unconditional 
stability.  The  solution  by  means  of  an  iterative  technique  is  the  main  topic 
of  this  paper;  moreover,  the  efficient  solution  by  means  of  the  Vale  sparse 
technique  will  be  the  main  result. 


DIFFERENCE  EQUATION  FORMULATION  OF  THE  3D  WIDE  ANGLE  PE 


We  are  concerned  with  the  solution  of  the  3D  wide  angle  PE,  Eq.  (4).  We 
seek  such  solution  by  means  of  the  Yale  sparse  technique,  in  particular, 
Craig's  method.  To  deal  with  the  solution  of  Eq.  (4),  we  must  first  discuss 
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the  solution  procedure  as  to  how  to  bring  Eq.  (4)  into  the  finite  difference 
equation  such  that  it  is  in  an  easy  and  acceptable  form  for  Craig's  method. 
Before  this  formulation,  we  have  a  few  definitions  to  state. 


Let  m  indicate  the  index  in  the  z  direction;  az  =  h  indicates  the 
z-increment.  Similarly,  is  used  to  indicate  the  index  in  the  s-direction; 
ao  =  s  is  the  ©-increment;  k  is  used  to  indicate  the  range  step  ar;  and  n  is 
used  to  indicate  the  range  level.  Also,  for  brevity,  define 


(6) 

(7) 

Then,  Eq.  (4)  can  be  expressed  in  a  short  expression  using  the  above 
definitions,  i.e.. 


x  =»  n2(r,a,z)  -1+4- 

ko  « 


and 


(kQr)  39 


Write 


3  /  i  *  Pi*  +  p?y\ 

—  u  =1  j-ik  +  ik  t— r — = — x — —  1  u 
ar  \  0  0  1  +  +  d2y  / 

_  /  l  +  Pi*  +  p?y\ 

■  po  *  ik0  rr-qp-f^yj  ■ 


then,  Eq.  (8)  can  be  written  in  a  short  operator  form,  i.e., 
IF  u  =^u  • 

Numerical  solution  to  Eq.  (10)  can  be  expressed  as: 

k_i 

n+1  3r  n 
u  =  e  u 


w 


(9) 


(10) 


(U) 
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Using  a  half-half  splitting  of  the  exponential,  and  setting  up  the  solution  to 
Eg.  (11)  uy  the  Crank-Nicolson  scheme  (an  implicit  finite  difference  scheme), 
we  find  that  an  implicit  finite  difference  discretization  to  Eq.  (4)  becomes 


[l  -  \  k*]un+1  =  [l  +  £  to*]un  . 


Using  the  definition  of*,  x,  and  y,  Eq.  (12)  becomes 


1  +  q  (n2(r,e,z)  -  1  +  V  +  q?  - J  un+1 

L  1  '  k2  3Z/  2  k2  (r+k )2  ae2J 

-  -  Ik  k  (p.-q.  )  /n2(r,a,z)  -  1  +  ^  +  (p?-q?)  -y — - *■  -Ap  u"+1 

20  1  1  V  k2  3z7  2  2  k2  (r+k  )2  ao" 

o  o 

=  1  +  di  ^n2(r ,o,z)  -  1  ♦  -y-  +  do  ~y~V  u0 

\  kg  az2/  2  k^  r  ae2  J 

(  U 

*  7  ik0  k  PrV  ("2(W)  -  1  4  7  75)  4  77?  ~2  u"  • 


2  2 

Using  central  differences  for  both  operators  — and  in  Eq.  (13),  and 

3Z2  302 


simplifying,  we  obtain 


J  +  qi(n"'1}  '  k2  h2  "  k2  (r^k)2  a2  '  2  k°  k(Prql)(n2_1) 


+  ^7  (Pl_ql)  +  \  (rik)2  a2 


4-11 


TD  7145 


*  (?  7  -  7  ^  ?)  *  (?  7  '  ^  lPl"qi ’hO 

(Prt>  777  7) 

*(|777  7 '  ^ (Prtl  777  7) 

.  (1  *  ,l(n2-D  -  ^  4  -  ^  4  4  *  7  ko  uprtX"2-1* 

V  1  »>  v  1 

lPl",l)  i?)1'1"*1-'  (Prai>7')u’-1-‘ 

t(|77‘H(P2'<'2)77)','"-w 
*(777*  7^ ‘^z1  77)  “«.*-1  ' 

This  is  the  large,  sparse  system  we  want  to  solve  efficiently. 


Let's  use  some 


abbreviated  symbols  to  simplify  the  coefficient  in  (14) 


Define 


/  2  2qll  2q 

V.-pi'"  -l)-7?7"7 


44(prV-^^7(Prt> 


2 
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n  ql  1  ,  1  k  ,  „  .  1 

^  jj-’fit  (prqi}  -j  * 


ko* 


1  \  1  k  /n  „  1  1 

^  O  L>  (  ) 


k£  (r+k)2  <s2  2  ko  22  (r+k)2  <s2 


■I.,-  (•  'W-H  -  ^7- -  ♦'[**. 

7^7  • 


k  1  /„  -  • 

”  k  ?  ' 

*o  ir  1 


_k 
o  r 


+  q2  1  1  .  .  1  k  /  _>  1  1 

R  =  -7  -7  -7  ♦  i  7  r—  (Pp-do)  ~7  ~7 

kf  r£  f  £  Ko  *  c  r  / 

o 


(15) 


We  can  see  Eq.  (15)  in  a  simpler  form,  i.e., 
P 


m,s,  in, a 
.+  n 


un.+1  +  QuJJ*}^  +  Qu"'’  -  +  R^'-i  +  Ru 


n+i 

m-l,a 


jn+1 

Jm,a+1 


n+1 

m,a-l 


=  W  V*  +  Q*  Cl, a  +  **  Cl, a  +  R+  C+l  +  R+  C'a-1 


(15) 


where  Q*  means  the  complex  conjugate  of  Q.  Pm  £  and  £  matrices 
depend  upon  the  variation  in  r,»,  and  z.  Q  matrix  is  constant  in  all  3 
variables.  R  and  R+  matrices  are  dependent  on  the  range  variable  only. 


AN  ILLUSTRATION 


For  illustrative  purposes,  we  use  a  simple  example  to  display  Eq.  (16)  in 
a  matrix  form.  In  general,  m  =>  1,  2,  ...,  M  and  =  1,  ...,  L.  Note  that 
m  =  0  indicates  the  surface  boundary;  m  =  M  +  1  indicates  the  bottom  boundary; 
and  8.  =  0  and  i  =  L  +  1  wi  1 1  be  explained  below. 
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We  start  with  assigning  m  =>  1,  2,  3,  4  and  a  -  1,  2,  3  at  the  initial 
range  level  n  and  march  to  the  next  range  n  +  1.  In  this  example,  L  =  3  and 
H  a  1,  2,  3  mean  there  are  three  sectors,  as  shown  in  Figure  1.  In 
computation,  we  must  deal  with  the  indexes  a  =»  0  and  i  a  4.  Since  the  index 
is  periodic  with  a  period  L  -  3,  then,  I  can  be  regarded  as  i  =  i  (mod  3). 
Therefore,  i  =*  0  is  the  same  as  t  -  3,  and  i  -  4  is  the  same  as  1,  as  also 
shown  in  Figure  1. 


Figure  1.  Azimuthal  Sectors 

Now  we  make  an  attempt  to  put  Eq.  (16)  in  a  matrix  form  encountering  the 
boundary  onditions.  We  use  the  convention  l  =  l  (mod  3)  to  express  n  and 

<4  by  £3  and  £3 ,  respectively.  We  then  can  construct  a  matrix  in  a 
general  form  making  use  of  the  periodic  boundary  condition.  We  find 
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In  general,  the  large  sparse  system  to  be  solved  is  in  the  form 
A  un+1  =  3  un  +  uf1  +  u"  , 

where  u"+1  contains  surface  and  bottom  boundary  information  at  the 


advanced  range 


level  and  u£  contains  surface  and  bottom  boundary 


information  at  the  present  range  level.  The  A  and  B  matrixes  possess  the 
format 

T  0  0  ...  0  0 

0  T  D  ...  0  0 

0  D  T  ...  0  0 

• 

0  0  0  .!.  T  D 

0  0  0  ...  0  T 

R  0  0  ...  0  D 

All  the  olock  matrices  (T,  0,  and  R)  are  of  the  same  order  MxM.  Each  T 

matrix  is  tridiagonal;  whereas  each  off-diagonal  block  matrix  D  and  R  are 

diagonal  matrices.  Entire  matrices  are  a  7-diagonal  matrix  with  the  property 
that  A  =»  At  and  B  »  BT. 

The  right-hand  side  of  Eq.  (18)  can  be  carried  out  by  one  matrix-vector 
operation  and  two  vector  additions.  Eq.  (18)  is  a  large,  sparse  system,  which 
we  want  to  solve  by  taking  advantage  of  the  Yale  sparse  technique. 

Note  that  if  we  consider  that  the  wave  propagates  all  around  a  complete 
360* ,  we  deal  with  a  system  where  A  and  B  are  of  the  form  (19),  i.e.,  a 

7-diagonal  matrix.  If  we  consider  that  the  wave  propagates  only  in  a  sector, 

then  the  periodic  boundary  condition  for  the  azimuthal  plane  disappears,  and 
we  then  solve  system  (18)  where  A  and  B  are  in  a  simpler  form  as  below: 
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T  0  0  ...  0  0  0 

D  T  D  ...  0  0  0 

0  0  T  ...  0  0  0 

• 

0  0  0  I  0 

0  0  0  ...  13  T 

0  0  0  ...  0  D 

which  is  a  5-diagonal  matrix. 


(20) 


Further,  if  we  consider  that  the  wave  propagates  only  in  a  vertical 
plane,  this  reduces  to  a  two-dimensional  case.  We  then  deal  with  the  system 
(18)  where  A  and  B  are  tridiagonal  matrices,  i.e.. 


T  0  0  ...  0  0  0 

0  T  0  ...  0  0  0 

0  0  T  ...  0  0  0 

0  0  0  ...  T  0 

0  0  0  ...  0  T 

0  0  0  ...  0  0 

It  is  important  to  note  that  when  p^  a  and  =  q2  =  0, 

the  system  (16)  reduces  to  the  3D  PE. 


0 

0 

T 


(21) 
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NUMERICAL  RESULTS 


As  a  test  of  accuracy,  the  Yale  sparse  technique  (Craig's  method)  was 
programmed  on  VAX  11/780  computer  to  solve  Eq.  (9)  using  the  system  of 
equations  expressed  by  Eq.  (16).  We  used  a  known  exact  solution  below  as  an 

accuracy  check. 


To  describe  the  test  procedure,  we  express  Eq.  (9)  in  the  form  of  Eq. 
( 5 ) ,  i  .e . , 


1  +  [(nztr.»,2)  -  1) 

O  u 


ik_ 


(p^) 


ar 


1 


(n2(r ,o,z )  -  l)  +^^]+  (p2-q2)7TT^  ,02  *  U 
0 


(k„r)  so 


We  look  for  a  solution  to  Eq.  (5)  in  the  form 


u(r,o,z)  =•  sin(fiz)  eimS  4>(r)  . 

Substituting  Eq.  (22)  into  Eq.  (5),  we  find 


(2 


(23) 


(p2-q2),n2  ! 
- - 7  '  U 
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We  select  n2(r,e,z)  -  1  -  %  =>  0  for  computational  simplicity.  Since 

ko 


<  =.  krtn(r,&,z)  =  w/c,  then 


*.  ■  tif  -  °f2  ■ 

Eq.  (23)  can  be  simplified  using  the  kQ  defined  by  Eq.  (24)  to  give 


li  f~ik°(P: 
8r  \  1  - 


-1UP  , 

2,n.  .. 


q2nf /(kQr) 


:) 


which  is  a  first  order  ordinary  differential  equation. 


The  solution  to  Eq.  (25)  can  readily  be  expressed  in  the  form 


iff(r)  dr 

$  =  A  e 


(25) 


(26 


The 


effort  needed  to  find  the  *(r)  is  the  evaluation  of  the  ff(r)  dr.  We  use 


u(rQl&,z)  =  sin(fiz)eirTI9  *(rQ) 

as  the  initial  field; 

u(r,o,zQ)  -  sin(flz0)eime  *(r)  =  0 

for  the  surface  condition  boundary;  and 

u(r,o,z8)  =  sin(nz8)eime  *{r)  -  0 


(27 

(28 

(29 
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for  the  bottom  boundary  condition.  These  boundary  conditions  are  particularly 
selected  such  that 

nz0  „  o  and  fizg  =  an  integer  multiple  of  it. 

The  initial  range  is  selected  to  start  at  50  m  so  that  the  farfield 
approximation  is  valid. 

The  azimuthal  plane  is  divided  into  10  sectors  at  36*  each.  Since  there 
are  10  sectors  and  we  partition  the  depth  into  199  increments,  then  we  solve  a 
system  of  equations  of  the  size  1990  x  1990  dealing  with  a  7-diagonal  matrix. 
The  results  presented  below  are  a  display  of  boundaries  between  two  adjacent 
sectors.  Not  only  do  we  compare  the  actual  computed  numerical  complex  numbers 
with  the  exact  solution  but  the  dB  values  as  well. 

Case  1:  Small  Angle  propagation  (p^  -  p2  -  1/2,  -  q2  »  0) 

An  evaluation  of  the  Jf(r)  dr  gives  -m2/(2k0r).  This  produces  the 
solution 


«J(r) 


2k  r 


A  e 


(31) 


Table  1  describes  the  results;  the  first  row  indicates  the  computed 
values  and  the  second  row  indicates  the  exact  solution.  The  results  are  taken 
at  the  boundary  between  the  third  and  the  fourth  sectors  at  108°  at  a  range  of 
50.4  in. 
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Table  1.  Results  of  Small  Angle  Propagation 


I 

...  z(D 

LOSS 

_  u(I) _  _ 

3 

30.00 

12.636 

(0.18834E+00 

-Q.13793E+00) 

3 

30.00 

12.636 

(0.18886E+00 

-0.13722E+00) 

6 

60.00 

6.859 

(0.36627E+00 

-0.26824E+00) 

6 

60.00 

6.859 

(0.36729E+00 

-0.26685E+00) 

9 

90.00 

3.749 

(Q.52397E+00 

-0.38372E+00) 

9 

90.00 

3.749 

l 0. 52541 E +00 

-0.38174E+00) 

12 

120.00 

1.841 

(0.65270E+00 

-0.47800E+00) 

12 

120.00 

1.841 

(0.65451E+00 

-0.47553E+00) 

15 

150.00 

0.688 

(0.74537E+00 

-0.54587E+00) 

15 

150.00 

0.688 

(0.74743E+00 

-0.54304E+00) 

18 

180.00 

0.108 

(0.79685E+00 

-0.58357E+00) 

18 

180.00 

0.108 

(0.799Q6E+00 

-0.58G55E+00) 

Case  2:  Wide  Angle  Propagation  (p^  =»  p2  *  3/4,  »  q2  »  1/4) 


An  evaluation  of  the 


dr  gives  the  solution 


<Hr) 


m(p2-q?)  k  r-m/gi 

-i  - -  an  - - 

2y€>  k  r+m£ 

A  e  c  0  c 


(32) 


Numerical  results  are  presented  in  the  same  manner  as  in  Case  t.  These 
results  are  given  in  Table  2. 
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0.117 

0.108 


G.89377E+00 

0.89492E+00 


(0.95544E+00 

(0.95673E+00 


-0.16245E+00) 

-0.16134E+00) 


-0.20236E+00) 

-0.20098E+00) 


-0.23150E+00) 

-0.22952E+00) 


-0.24602E+00) 

-0.24537E+00) 


From  the  solution  results,  we  examine  the  behavior  of  the  solution  of 
Case  2  for  large  kQr.  First,  we  consider  the  real  part  of  the  solution  for 


in(P?-q?)  fk  r-mvC 

x  „  — 5 — £~  anl  —S; - - 


rv^l  i 

Lk0r+inrT2J 


m(Pp-qo)  i  nrq„ 

cos(x)  =  cos - — —  --XT"  7 - 1 

2^2  V  k  (k^ 


m^2  1  m  q2 
V  “  ?  0cor)2  *** 


cos(x)  s  cos 


m(p?-q?)  T  2 /Join'll  nr(p„-q?) 

- — —  -  T7-  =  COS - r— - 
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for  large  kQr.  For  wide  angle  parameters  (33)  reduces 

to  the  real  part  of  the  solution  of  Case  1. 

In  examining  the  cos  ■?)>  we  note  that  the  function  increases 

2  ® 

monotonical ly  after  kQ  «  m  /*.  It  approaches  unity  as  kQr  •>  <*>.  When 
the  function  is  close  to  unity,  three-dimensional  effects  are  lost,  and  Eq. 

(4)  behaves  like  the  two-dimensional  parabolic  wave  equation  below 


2 

The  more  rapid  the  azimuthal  variation  (i.e.,  the  bigger  m  ),  the  further 
out  in  range  the  three-dimensional  effects  influence  the  solution. 

An  application  is  also  presented  here  as  a  second  example.  This  example 
has  peen  solved  by  Perkins  and  Baer,11  using  the  Split-step  algorithm  for 
three  dimensions.  The  sound  speed  profile  is  taken  from  a  Pacific  profile 
such  that  c(r,e,z)  «  cm(z)  +  (O.OOl)rsine,  where  cm(z)  takes  on  the 
tabulated  values  below. 
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Table  3.  A  Pacific  Profile 


z(m) 

c(z)  (m/s) 

0.00 

1536.500 

152.400 

1539.243 

406.300 

1501.143 

1015.9 

1471.882 

5587.91 

1549.606 

5587.91 

1555.526 

This  profile  has  a  large  linear  gradient  in  the  cross  range  direction; 
the  gradient  is  1  m/s  per  km.  The  profile  in  the  vertical  plane  at  0*  is  a 
typical  profile  in  the  North  Pacific  Ocean. 


The  source  is  placed  at  254  in  below  the  surface  with  a  source  frequency 
of  25  Hz.  We  calculated  the  propagation  loss  up  to  a  maximum  range  of  140  km. 
We  choose  to  present  below  the  results  on  one  particular  sector  at  0*.  Along 
with  the  30  wide  angle  PE  solution  plot  is  the  graphical  result  of  Perkins  and 
Baer^  for  comparison.  The  propagation  loss  reading  at  120  km  for  the  same 
receiver  depth  is  approximately  90  dB,  showing  satisfactory  agreement  with  the 
known  result.  The  30  wide  angle  PE  result  is  presented  in  figure  2;  and  the 
Perkins-Baer  result  is  presented  in  figure  3. 
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CONCLUSIONS 


In  situations  where  FFT  is  applicable,  the  FFT  can  do  well.  Under  the 
same  situation,  the  CG  method  can  also  solve  the  same  problem  with  the  same 
accuracy,  but  the  computation  speed  is  not  competitive  with  the  FFT 
computation.  The  solution  to  the  wide  angle  three-dimensional  partial 
differential  equation  cannot  be  directly  solved  by  the  FFT;  this  is  a  definite 
advantage  of  the  Yale  sparse  technique.  The  applicable  CGNE  that  we  used  here 
requires  5N  multiplications  per  loop  plus  two  matrix-vector  products,  and  only 
4N  storage  locations  are  required  for  the  vector  operations. 

Since  CG  is  an  iterative  technique  dealing  with  inner  products,  it  is 
desirable  to  implement  the  procedure  in  a  vectorized  machine.  This  is  another 
advantage  of  the  Yale  sparse  technique. 

The  numerical  solutions  produced  in  this  paper  demonstrated  the  general 
purpose  capability  of  the  Yale  sparse  technique.  Even  though  the  solution  is 
accurate,  the  present  solution  can  by  no  means  be  regarded  as  the  most 
efficient  solution.  It  is  believed  that  a  clever  preconditioning  technique 
can  be  developed  to  enhance  the  efficiency  of  our  applications. 
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5.  A  RANGE  REFRACTION  PARABOLIC  EQUATION 


Frederick  D.  Tappert 
University  of  Miami 

Ding  Lee 

Naval  Underwater  Systems  Center 


ABSTRACT:  Application  of  the  standard  parabolic  wave  equation  to 
solve  real  problems  requires  a  clever  selection  of  the  reference 
wavenumber  k0.  An  extended  parabolic  equation  (PE)  having  range 
refraction  capability  Is  reintroduced  to  be  totally  Independent  of 
k0.  The  existing  implicit  finite  difference  (IFO)  model  was 
applied  to  test  the  range  refraction  PE.  Results  compare  favorably 
with  known  solutions  for  weakly  range-dependent  environments,  but 
yield  significant  corrections  for  propagation  through  strong 
oceanic  fronts. 


5-1/5-2 
Reverse  Blank 
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INTRODUCTION 


The  RAnge  REfraction  Parabolic  Equation  (RAREPE),  introduced  by  Tappert 
[1]  over  a  decade  ago,  has  important  acoustic  effects  but  has  been  neglected. 
Since  the  standard  parabolic  equation  (PE)  is  performing  satisfactorily,  users 
have  not  paid  attention  to  the  RAREPE.  Besides,  there  did  not  exist  efficient 
algoritnms  directly  applicable  to  solve  the  RAREPE.  Now,  we  have  the  implicit 
finite-difference  (IFD  [2])  package,  and  the  effort  required  to  modify  the 
available  IFD  code  to  solve  the  RAREPE  is  inexpensive.  To  give  a  complete 
understanding  of  the  RAREPE,  we  first  summarize  the  derivation  of  the  RAREPE, 
then  describe  how  we  solved  the  RAREPE  by  the  finite  difference  solution.  A 
special  section  is  devoted  to  discuss  a  set  of  illustrative  examples.  These 
examples  are  used  to  show  (1)  the  close  agreement  between  the  standard  PE  and 
the  RAREPE  if  there  is  no  front,  (2)  the  important  property,  independence  of 
kQ,  of  the  RAREPE;  and  (3)  the  range  refraction  effects  by  weak,  moderate, 
and  strong  fronts. 


DERIVATION  OF  THE  RANGE  REFRACTION  PE 

We  start  with  the  two-dimensional  reduced  wave  equation,  i.e., 

Prr  +  y  Pr  +  Pzz  +  k l  n2(r,z)p  .  0  (1) 

Setting  p(r,z)  =  ’-z ^  and  applying  the  farfield  approximation,  kQr  »  1,  we 

find  Eq.  (1)  becomes 

V  +  “zz  +  kon2(r.*)u  =  0  .  (2) 

Write  Eq.  (2)  in  the  form 
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&  *  k“q2)“' =  ° 


where 


q2  =  n2(r,z)  +  . 

k2  3Z2 


(3) 


(4) 


We  assume  that  the  dependence  of  n(r,z)  on  the  range  variable  r  is  weak 

such  that  an^-z^-  is  negligible,  but  we  shall  later  pick  up  the  neglected 

a  2 

terms.  This  assumption  allows  the  operator  —  to  commute  with  Q  .  We  can  now 

9  i 

factor  Eq.  (3)  into  two  equations: 

3U. 


and 


1  ~  +  ko  QQ+  "  0  * 


3U 

-1  ST  +  ko  Qfi.  “  0  • 


(5) 


(6) 


Then,  the  solution  field  u(r,z)  is  just  the  combination  of  the  outgoing  wave 
u+(r,z)  and  incoming  wave  u_(r,z). 

The  operator  q,  defined  by  Eq.  (4),  is  actually 


Q  -  (1  +  (n  (r,z)  -  1)  +  -*• 

kz 


1  a2  V 

IT?) 


12 


(7) 


In  this  paper,  we  deal  with  the  small  angle  PE,  i.e.,  we  approximate  the  Q  by 


Q  =  1 


(n2(r,z)  -  1)  +  — \  o 
k2  3ZZ 


(8) 


ikQr 


and  make  use  of  the  usual  envelope  definition  u  =  u(rvz)  e  ,  which  leads  to 


5-4 


TD  7145 


the  standard  PE  originally  introduced  by  Tappert  [1]»  i.e- 


To  make  the  local  error  small,  we  need 


|n2(r,z)  -  1| |  «  1  , 


Ik2  az2l 
0 


2  1  II 

A  detailed  discussion  of  the  estimate  of  ||n  (r,z)  -  1||  and  ||-s-  jll  can  be 

k0  32 

found  in  reference  1. 


To  keep  track  of  the  relative  errors  made  in  the  course  of  calculation, 

1  *2 

we  can  monitor  the  size  of  ||n2 (r,z)  -  1 1|  and  ||-* — -yll  and  keep  them  both 

ko  32 

small.  We  now  present  a  modified  PE,  which  requires  that  only  one  of  the 
n2(r,z)  -  1  and  be  small,  as  occurring  in  Q,  but  is  of  order  unity. 

lko  32 

This  modified  PE  is  capable  of  dealing  with  a  large  range  variation  of  the 
index  of  refraction. 


2  1  ^ 

Consider  two  operators  A  and  B, where  A  =  n  (r,z)  and  B  =  -j — t.  In 

k  3  z 
0 

1/2 

general  A  and  8  do  not  commute.  We  expand  (A  +  60)  by  the  formula 
(A  +  6B)1/2  =  A1/2  +  «C  +  o(62)  , 
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where 


co 


o 


The  proof  of  (10)  is  given  in  reference  1. 

We  now  apply  Eqs.  (10)  and  (11)  to  the  operator  Q  given  by  (7). 
1  a2 

For  small  — r  — y,  we  obtain 


It  can  oe  found  that 


00 


Substituting  Eq.  (13)  into  Eq.  (5)  for  the  outgoing  wave,  we  obtain 


Eq.  (14)  is  the  RAREPE  and  is  valid  to  all  orders  in  (n2(r,z)  -  1). 

Equation  (14)  is  approximately  equal  to  the  standard  PE,  Eq.  (9),  when  n(r,z) 
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is  range-independent  and  is  close  to  unity;  in  this  case  the 
in  the  standard  PE  can  be  replaced  by  (n(r,z)  -  1). 


n^(r.z)  -  1 


It  ought  to  be  noted  that  Eq.  (14)  and  Eq.  (9)  are  the  same  if  the  n(r,z) 
is  a  constant  in  both  r  and  z  variables,  and  n(r,z)  =  1.  Substituting  the 
constant  n(r,z)  into  Eq.  (14),  we  find 


2 

3U  ,  >  .  1  3  U 

,k0  n(r,2)  U  ♦  5i^ 


Eq.  (15)  can  be  put  in  the  form 


and  replacing  (n(r,z)  -  1)  by  (l/2)(n2(r,z)  -  1),  we  find 


$  -  7  k0(n2(r,z)  -  D^+ 

0  3  Z 


We  see  that  Eq.  (17)  is  exactly  in  the  same  format  as  Eq.  (9). 


We  now  have  established  the  relationship  between  (15)  and  (16)  based  on 
the  transformation 


u(r,z)  =  i|»(r,z)  e  .  (1* 

We  now  proceed  to  show  that  solutions  to  Eq.  (15)  and  Eq.  (16)  are  identical 
in  magnitude  in  order  to  establish  the  k  independence. 
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Substituting  Eq.  (18)  into  Eq.  (15),  we  find 

-  ikQ(n(r,z)  -  1)*  F  *zz  ’  (19) 

which  is  identical  to  Eq.  (16).  From  the  relationship  of  Eq.  (18),  it  is 
easily  seen  that 


This  shows  that  the  solution  of  (15)  is  equivalent  to  the  solution  of  (18)  in 
magnitude  independent  of  kQ. 

THE  REFERENCE  WAVENUMBER  kQ 

We  return  now  to  the  standard,  small  angle  PE,  which  takes  the  form 

2  k 

1  ar-  +  W~  "Hr  +  ~7  (n2(r»z)  -  Du  -  0  .  (20) 

0  3Z 

It  is  clear  that  Eq.  (20)  is  kQ-dependent,  thus,  different  kQ's  lead  to 
different  solutions.  Obviously,  there  is  only  one  kQ  associated  with  a 
given  set  of  environmental  conditions  that  will  produce  the  solution  closest 
to  the  real  solution.  We  have  been  confronting  the  problem  of  how  to  select 
the  best  kQ.  in  fact,  PE  users  never  have  to  worry  about  the  kQ-selection 
because  existing  PE  models,  such  as  the  split-step  code  [3]  and  the  IFD 
package  [2],  all  offer  the  option  to  have  a  default  kQ  value  if  the  user  is 
not  certain  what  kQ  to  use.  The  default  kQ  is  chosen  from  kQ  = 
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2’rf/Co,  where  c0  is  the  reference  speed  and  is  selected  to  be  the  average 
sound  speed  in  the  water  column.  The  user  can  choose  kQ  to  be  the  sound 
speed  at  the  source  level  for  range-independent  environments.  Moreover,  for 
range-dependent  environments  and  for  range-dependent  sound  speed  profiles,  the 
user  can  ask  that  an  interpolation  of  the  sound  speed  profile  be  performed 
within  each  range  interval  and  apply  the  same  procedure  to  select  cQ  as  in 
the  range-independent  case.  These  choices  so  far  present  no  big  problem;  and 
make  the  selection  of  the  kQ  ignorable.  Pierce  [4]  re-emphasized  the 
importance  of  the  kQ  selection  and  introduced  a  formula  to  determine  the 
range-dependent  kQ  based  on  the  Rayleigh  quotient.  Some  numerical 
experiments  have  been  carried  out  at  NUSC,  New  London  Laboratory.  The  results 
show  some  phase  effects  on  kQ  variation.  A  detailed  study  of  the 
kQ-selection  is  going  to  be  reported  separately  when  it  is  completed.  All 
these  facts  strongly  suggest  the  desirability  of  having  either  a  variable  kQ 
PE  selection  or  a  kQ-independent  PE.  This  paper  chooses  to  deal  with  the 
latter. 


THE  FINITE  DIFFERENCE  FORMULATION  AND  SOLUTION 


Rewrite  Eq.  (14)  in  the  form 


i  +  _L  _L  (l  _l\  +  k 

ar  2kQ  az  az  J  o 


vu  =  0 


where 


v  =  n  + 


1  (  nz2  nzz 

4k  1  ■“ J  "  T 
o  \n  n 


(21) 


(22) 
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vm  '  f  *  4k2  n4  7  '  Vl  Vl)  ' 

^  *0  m  -* 

c 

Using  tne  definition  k  n  =  7 - 7  =  7=  k,  (27)  can  be  written  as 

0  cQ  c  c 

'  *  ^  y  '  vTvr)  ‘ 

Substituting  (24)  through  (28)  into  Eq.  (21),  we  find 


(27) 


(28) 


30, 


m 


ar 


=  -  at  U 


U.,+B  U  +  Y  U  , 

m  m+1  m  m  Tm  m-1 


(29) 


where 


fC 

Note  that  a  =  -a „  and  v  =  -ym  , .  We  see  that 
m  in  m  m-1 

=  am(um  um+l  -  um+l  um)  +  Ym(um  um-l  _  um-l  um) 


3 1..  i2 
ar 


This  implies  that 


~(T\U  |2>)  =  V  a  (u*  u  „  -  u*  ,  u)  -  V\  (u*  u  .  -  u*  .  u  ) 

arl/  jl  m1  J  /  1  m\  m  m+1  m-1  m/  /  ^'my  m  m-1  m-1  m/ 

x  m  '  m  m 


Preceding  Page  Blank 
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C°m  ~  Ym+l)(um  um+l  "*  ^m+l  um) 


But’  “m  =  Ym+l*  Therefore 


a  r 


(EwV°>  Ei 


uJ  a  constant 
ml 


This  implies  that  the  finite  difference  scheme  is  energy  conservative. 

Next,  from  Eq.  (29),  we  set  up  the  Crank-Nicolson  difference  equation  as 
follows: 


k  n+1  n+1  .  /.  k  0n+l\  n+1  k  n+1  ,.n+l 

7  um+l  +  \X  -  7  Bm  )  um  "  7  Ym  Vl 

k  n  ,,n  +  /,  x  k  „n\  ,,n  ,  k  n  n 

7  “in  Vl  +  \l  +  7  Bmj  um  +  7  Ym  um-l  * 


(30) 


This  is  the  exact  IFD  format  recognized  by  the  IFD  model.  The  solution  by  the 
IFD  code  becomes  easy. 


NUMERICAL  ILLUSTRATIONS 

For  qualitative  information,  this  section  presents  three  examples  which 
are  used  to  show  the  various  effects  of  range  refraction. 

We  use  the  following  cannonical  profile  whose  input  parameters  are 
defined  and  tabulated  below: 
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c(r,z)  =  cA[l  +  e(e~n  -  1  +  n)J 

tfiere 

2  -  ZA 

n  =  (8/2 )  * 

=  sound  speed  of  axis, 

z  =  depth  variable, 

z A  =  depth  of  axis  of  sound  channel,  and 
B  =  thickness  of  thermal  front. 
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(32) 


We  assign  an  ocean  depth  of  5  km  and  assume  the  ocean  bottom  to  be  flat. 
We  calculate  the  propagation  loss  up  to  140  km  in  range.  We  place  the  source 
at  100  m  below  the  surface  with  a  source  frequency  of  100  Hz. 


Define 


B(r)  =  Bx 


where  r^  is  the  range  at  which  the  front  occurs, 
L  is  the  length  of  the  front,  and 
Bj  are  parameters. 


In  addition,  define 


(33 


e(r) 


(34 


where  g  =  2cA  /B 


(35 
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PROFILE 

B1 

(km) 

82 

(km) 

(km) 

rF 

(km) 

(km) 

REMARKS 

1 

1.2 

1.2 

1 

50,  60 

_ 

No  front 

2 

1.2 

1.0 

1 

50,  60 

20 

moderate 

front 

3 

1.2 

0.8 

1 

50,  60 

20 

strong 

front 

The  cannonical  profiles  for  these  different  examples  are  described  in  figure  1. 

SOUND  SPEED  (m/s) 


Figure  1.  Cannonical  Profiles 

The  examples  (using  the  above  three  profits)  are  executed  using  the  following 
information: 
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EXAMPLE 

INITIAL 

RANGE 

(kin) 

INITIAL 

PROFILE 

1 

1  1 

FINAL 

PROFILE 

RANGE 

FRONT 

(kin) 

1 

0 

1 

140 

1 

No 

2 

0 

1 

140 

2 

60 

3 

0 

1 

140 

3 

60 

The  following  set  of  graphs  (figure  2)  were  obtained  by  the  IFD  model 


using  the  input  information  of  Profile  1. 


lots 
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The  results  have  the  following  meanings: 


CASE 

RANGE  FRONT  AT 
km  RANGE 

FRONT 

SOLVED  BY 

A 

No 

Standard  PE 

3 

60 

Moderate 

RAREPE 

C 

60 

Strong 

RAREPE 

D 

50 

Strong 

RAREPE 

Figure  3  presents  the  propagation  loss  curves  over  the  range  interval  [0, 
140  km]  for  three  different  receiver  depths.  There  is  no  front  present  in 
these  examples.  The  left  column  displays  the  standard  PE  results;  the  right 
column  displays  the  results  with  range  refraction.  Notice  that  the 
differences  among  these  results  are  very  small.  In  order  to  make  the 
difference  more  visible,  we  must  group  the  results  together  in  a  magnified 
plot. 
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0  20  40  60  80  100  120  140  0  20  40  60  80  100  120  140 


0  20  40  60  80  100  120  140  0  20  40  60  80  100  120  140 


RANGE  (Krr;. 

Figure  3.  Results  Without  Range  Front 

The  original  Helmholtz  equation  is  ^-independent.  For  a 
range-inuependent  problem,  we  use  an  accurate  fast  field  program  iFFP  [3]) 
solution  as  a  reference  for  comparison.  The  split-step  [1]  solution  is  also 
included.  Figure  4  contains  the  solutions  produced  by  the  1FD,  the 
split-step,  aria  the  FFP,  with  ana  without  range  refraction  effects. 


Figure  4.  Results  Comparison 
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From  the  comparison,  a  little  difference  is  shown  among  the  IFD, 
split-step,  and  the  FFP  solutions.  When  comparing  results  among  IFD,  IFD 
range  refraction,  and  the  FFP,  also  a  little  difference  is  shown.  However, 
the  comparisons  among  the  IFD  wide  angle,  IFD  range  refraction,  and  the  FFP, 
we  experience  a  difference  between  the  IFD  wide  angle  and  the  IFD  range 
refraction;  this  is  expected  because  the  present  IFD  wide  angle  model 
accommodates  the  range  refraction.  This  is  going  to  be  pointed  out  in  the 
next  example. 

Similar  as  the  set  of  no  range  refraction  results,  the  next  set  (figure 
5)  consists  of  propagation  loss  curves  over  the  range  interval  [0,  140  km] 
with  a  strong  range  front  that  occurs  at  the  range  of  60  km. 
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Figure  5.  Results  With  a  Strong  Range  Front 

Whenever  trie  range  front  is  present,  especially  the  strong  front,  the 
RAREPE  results  are  more  meaningful  than  the  standard  PE  results.  Results  on 
figure  6  show  a  difference  between  the  RAREPE  and  the  wide  angle  PE.  Both  the 
standaru  PE  and  the  Wide  Angle  PE  were  formulated  dependent  upon  the  special 
approximation  of  the  square  root  operator ^1  +  e  +  p,  where  e  =  n 2(r,z)  -  1, 

1  3  2 

ana  p  =  —*■ — .  The  IFD  solution  considers  the  square  root  to  be  approximated 
ko  3Z 

in  the  form  { qj'e~  I  •  When  p  =  ^  dnu  q  =  0,  the  resulting  PE  is  the 

stanoaru  PE.  When  p  =  3/4  and  q  =  1/4,  the  resulting  PE  is  the  wide  angle  PE. 

Ex  panning  (1  +  q(e  +  y/'1  U  +  PU  +  y>)  =  (1  -  qU  +  p/  +  q2U  +  uj2 

+  •••/  (1  +  pU  +  UJ)  =  1  +  (p  -  qMe  +  u)  -  q IP  -  q ) ( e  +  u)2  +  high  order 

terms.  Notice  that  trie  portion  (e  +  y)2  represents  the  wide  angle  ana 
(e  +  y)2  =  e2  +  y2  +  ep  +  ye ,  where  the  last  two  terms  represent  the  range 
refraction.  Therefore,  the  present  wioe  angle  PE  IFD  version  accommodates  the 
range  refraction. 
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- IFD  Wide  ANGLE 

— - IPO  Range  Refraction 

. IFD  No  Range  Refraction 


Figure  6.  Results  Comparison 
CONCLUSIONS 


A  parabolic  equation,  having  range  refraction  capability,  is 
re-i ntroducea.  This  range  refraction  PE  is  independent  of  the  choice  of  kQ, 
ana  is  also  more  efficient  to  handle  a  rapialy  varying  index  of  refraction  on 
the  range  variable.  Numerical  results  show  that  the  range  refraction  PE  is 
useful  for  strong  range  fronts.  Under  such  environments,  the  range  refraction 
PE  has  an  aavantage  over  trie  standard  PE.  However,  a  price  has  to  be  paid  on 
tiie  computation  speed. 


b-k2 
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6.  THE  HYBRID  PARABOLIC  EQUATION  —  A  RAY  MODEL 
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ABSTRACT:  Using  the  standard  parabolic  equation  to  solve  high 
frequency  problems  Is  Impractical  because  of  excessive  running 
time.  A  new  HYbrld  Parabolic  Equation  using  Ray  theory  (HYPER)  is 
developed  to  handle  high  frequency  problems.  This  paper  discusses 
the  theory  and  development  of  the  HYPER  discussed  In  detail. 
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INTRODUCTION 

Successful  applications  of  the  standard  parabolic  equation  (PE)  [1]  have 
been  evident.  However,  using  the  standard  PE  to  solve  high  frequency  problems 
suffers  an  excessive  running  time  [2].  To  overcome  this  difficulty,  a  new 
HYbrid  Parabolic  Equation  using  the  Ray  theory  (HYPER)  is  designed  to  be 
particularly  effective  in  handling  high  frequency  problems.  A  complete 
discussion  is  given  to  describe  how  the  HYPER  is  derived.  Within  the 
discussion,  how  the  ray  equation  is  obtained  will  also  be  described  in  full. 

THEORETICAL  BACKGROUND  AND  DEVELOPMENT 

Prior  to  the  development  of  the  high  frequency  parabolic  equation 
(HYPER),  it  is  necessary  to  list  a  set  of  definitions  below. 

Let  u(r,z)  be  the  pressure  field,  n(r,z)  be  the  index  of  refraction,  kQ 
be  the  reference  wavenumber,  cQ  be  the  reference  sound  speed,  and  c(r,z)  be 
the  sound  speed  profile. 

The  standard  PE  takes  the  general  expression 


1  -  +  2T  *~T  “  ko  U(r’z)  u  =  0 

*0  3ZC 


ar 


(1 


for  small  propagation  angles,  and  U(r 


,z)  =  i(l-n2(r,z))  =  l  (l  - 
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An  observation  of  Eq.  (1)  stirs  up  two  motivations  that  lead  to  the 
development  of  the  HYPER 

Motivation  1:  The  standard  PE,  Eq.  (1),  has  no  limit  on  frequency,  but  the 

equation  itself  is  strongly  frequency  dependent.  As  a  function 
of  frequency,  the  computation  time  is  excessive;  because  of 
this  we  need  a  more  effective  way  to  handle  high  frequency 
problems.  We  feel  that  there  is  no  reason  why  we  cannot 
develop  a  high  frequency  PE  that  possesses  the  same  format  as 
(i). 

Motivation  2:  The  geometry  of  the  acoustics  is  independent  of  frequency.  It 
is  highly  desirable  to  have  a  ray  PE  model  independent  of 
frequency  but  with  full  range  effects.  We,  thus,  seek  to 
develop  a  high  frequency  PE  in  format  (1)  but  totally  frequency 
independent. 


6-4 


z  =  aR(r)  +  c  , 
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the  dot  indicates  the  r-derivative;  SQ(r)  is  taken  for  granted  as  a 
function  of  r  for  the  time  being  and  is  defined  later  in  the  section  on  the 
development  of  ray  equation  (the  particular  ray  equation  associated  with 
(2)).  Now,  we  proceed  to  show  how  (2)  was  derived.  Write 

ik  S(r,z) 

u(r,z)  =■  A(r,z)  e  .  (5) 


For  large  kQ,  substitute  (5)  into  the  standard  PE,  Eq.  (1),  and  use 
asymptotic  expansions  for  kQ  (keeping  order  up  to  o(l/k0)),  we  find 


(6) 


Equation  (6)  is  an  inhomogeneous,  nonlinear  partial  differential  equation 
(PDE)  called  tne  Hamilton-Jacoby  equation.  Rays  of  the  PDE  (6)  are  the 
characteristics  of  the  PDE  (6).  To  solve  eq.  (6),  differentiate  (6)  with 
respect  to  z,  we  find 


a2S  +  aS  a2S  +  aU 
araz  az  az 


0 


Then,  let  e 


Write 


but. 


aS/az  and  substitute  it  into  the  above  equation,  we  obtain 


a®  + 

e  2® 

+■ 

aU 

ar 

az 

3Z 

de 

89  . 

e 

39 

Hr  a 

ar 

3  Z 

de  _ 

38  * 

dz 

3fe 

dr 

ar 

Hr  az 

(7) 
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Comparing  the  above  2  equations,  &  must  equal  dz/dr.  Solving  (7)  by  the 
method  of  characteristics,  we  have  on  the  curve  z(r) 

dz 

HF  "  ®  » 

,  de  dU  n 

and  dr  dz*0  * 


Combining  (8)  and  (9)  gives  the  HYPER  ray  equation,  i.e.. 


-^U(r.z)  . 


Suppose  we  trace  a  particular  ray  (i.e.,  a  particular  solution  z  ■  zR(r)). 

To  deal  with  the  standard  PE,  let  5*  z  -  zR(r).  Corresponding  to  the 

operator  jp,  we  find  -  zR  ;  similarly  ~  for  —  .  The  standard  PE  can  be 
transformed  into 

1  fr  "  iiR  +  2T-  M  “  ko  U(r>zR(r)  +  t)  u  *  0  •  (11) 

0  35 

Now,  recall  the  SQ(r)-.  We  express 

S Q(r)  -  S(r,zR(r) )  ,  (12) 


and  we  refresh  our  memory  that  the  S(r,zR(r))  satisfies  the  Hami lton-Jacoby 
equation  (6). 

From  the  relationship  of  9,  we  see  that 

4  dS  dz 
*  3  Hz  3  W  • 
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S(r,z) 


■  /ffd2  */^dr  -/(ft  #  *  w) dr 


which  implies  that 


This  is  the  complete  expression  for  SQ(r).  We  need  the  r-derivative  of 
SQ(r).  We  find  that 

3FV>  -?(£*)  -  U(r.zR(r)} 


(13 


(14 


Now,  substitute  expression  (4)  into  Eq.  (2),  we  obtain  the  desired  HYPER,  i.e. 

-  ko  [u<r,*RCr)  *  c)  -  U(r,t|,(r» 

-tjf  U(r,zR(r))]  u  •  0  .  (1! 

Next,  we  shall  show  how  (15)  was  obtained  making  use  of  all  the  developments. 
From  (4),  we  find 
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a  li 
ar 


and 


au 

ac 


a2u 

17 


au 

ar 


g  +  ikft  2 
ac  o 


iko(?  !S  -  UR  *”«)“ 

R  ;]e,ko[so(r)t  tzR(r)] 

=] 


,1ko[Vr)*  'Vr)] 


2-  au  2  2  ~1  1ko[So(r)+  ?2R(r)] 

_*2  'koJRlF-ko<!R)2“h 


dy  substituting  the  above  partial  derivatives  into  (2),  we  obtain  (11). 


Q.E.D. 


COMPUTATIONAL  ASPECTS 


To  perforin  the  computation  completely,  a  number  of  steps  are  involved. 
Each  present  step  depends  upon  the  previous  step  after  the  problem  is  properly 
started.  We  discuss  the  computations  required  for  each  step  and  show  the 
continuity  from  one  step  to  the  next.  Most  of  the  present  computations  are 
carried  in  a  practical  manner  to  make  the  solution  work.  Room  exists  for 
future  improvement  in  computations. 


Step  1:  Calculation  of  the  Ray  Equation 

One  important  portion  of  the  computation  is  the  calculation  of  the  ray 
equation  to  describe  how  the  ray  is  traced.  The  ray  trace  portion  requires 
the  implementation  of  the  ray  equation,  (10),  i.e., 


6-8 


TD  7145 


aU 

3Z  ' 


The  ray  equation  is  a  second  order  ordinary  differential  equation  (ODE)  and 
can  be  treated  as  a  purely  initial  value  problem.  This  second  order  ODE  can 
be  solved  efficiently  by  an  existing  convergent  numerical  ODE,  equation  such 
as  the  methods  given  in  reference  3.  We  apply  the  Stomer-Cowell  formula  to 
perform  the  ray  trace. 


We  express  our  ray  equation  in  a  short  form  as 
z"  • 

The  Cowell  corrector  formula  takes  the  form 

-  22n  *  Vl  -  ^  [Vl  *  10fn  *  Vl]  • 

To  get  the  predicted  value  to  integrate  formula  (17),  we  use  the  Stormer 
predictor  formula 


P 

Vl 


Zz, 


Vl 


Ur)2  fr 


(16) 


(17) 


(18) 


Since  we  are  solving  an  initial  value  problem,  the  z(0),  z'(0)  are  known.  We 

also  know  fR  and  zq,  which  is  z(0).  We  need  to  know  zj.,  which  we  choose 

to  obtain  by  the  Taylor  series  expansion,  i.e., 

7.j  -  z (0)  +  Ur)  z'(0)  +  ^Ur)2  z"(0)  +  ... 

-  z  (0)  +  Ur)  z'  (0)  +  -j(ar)2  f  (0)  +  ... 
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The  Taylor  expansion  has  an  order  of  accuracy  in  truncation  error  less 
accurate  than  the  numerical  ODE  methods  used.  We  overcome  this  difficulty  by 
using  a  very  small  step  size. 

According  to  Henrici  [3],  the  selection  of  the  step  size  to  satisfy  the 

1  1/2 

corrector's  convergence  is  such  that  Ar  <  — j -  ,  where  L  is  the  Lipschitz 

I?  L 

constant  of  f(r,z). 

To  carry  out  the  complete  procedure  will  give  zR(r). 

Step  2:  Determination  of  the  Ray  Trace  Region 

Next,  we  make  use  of  the  information  obtained  from  the  previous  step,  and 
attempt  to  define  a  strip  (called  the  wide  W)  covering  the  ray  path.  The 
figure  below  shows  the  picture  in  the  r-z  plane. 
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here  zs  is  the  source  depth  and,  zr  is  the  receiver  depth.  The  region  of 
interest  is  determined  to  be 


zs-wiz£zs+w  * 

therefore,  the  width  is  2W.  We  see  that  if  we  solve  the  problem  in  the 
original  r-z  plane,  we  need  to  solve  the  problem  in  the  entire  original 
region;  therefore,  the  amount  of  work  is  by  no  means  able  to  be  reduced. 

The  advantage  is  evident  if  we  solve  the  problem  in  the  equivalent  r- 
plane  as  shown  by  the  figure  below. 


Zs 


r 


€ 


The  region  covered  in  the  r-£  plane  is  small;  therefore,  much  less  work  is 
needed  to  complete  the  computation  and  should  be  more  accurate.  This  above 
illustration  shows  one  single  ray;  for  a  family  of  rays,  we  handle  the  family 
by  taking  the  union. 
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Step  3:  Computation  of  V(r, t) 

Function  V(r,C)  is  calculated  using  Formula  (3),  i.e., 
V(r.c)  -  U(r,zR(r)  +  ?)  -  U(r,zR(r)) 

where  U(r,z)  is  defined  to  be 


On  the  ray  path,  ll(r,z)  is  calculated  by 


and  j|  U  is  calculated  by 

jiU(r,,R(r)).ji  [u(r,zK(r)  *  az)  -  U(r,zR<r))]  . 

Step  4:  Computation  of  SQ(r) 

SQ(r)  =.  S(r,ZR(r))  is  calculated  by  formula  (13). 

Step  5:  To  Obtain  the  Solution  u(r,z) 

Our  main  objective  is  to  obtain  the  solution  of  Eq.  (1),  i.e.,  u(r,z). 
We  summarize  the  procedures  involved  in  order  to  obtain  the  u(r,z)  and  show 
how  the  u(r,z)  is  obtained. 
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Obtain  the  ray  equation  solution  zR(r). 

Using  zR(r),  compute  V(r,?)  by  means  of  V(r,t)  -  U(r,z) 
relationship. 

Group  a  family  of  rays  together  that  have  strong  effects  between  the 
source  and  receiver. 

Determine  the  grid  points  in  the  r-c  plane  and  set  up  the  numerical 
solution  for  u(r,c)* 

Solve  u(r,t)  by  the  IFD  model. 

Solve  SQ(r)  »  S(r,zR(r))  by  formula  (13). 

Final ly. 


u(r,z) 


u(r,c) 


1ko[So(r)+  2R(r)] 
e 


CONCLUSIONS 


We  have  developed  the  HYPER  specially  for  handling  high  frequency 
problems.  Through  an  efficient  implementation,  the  advantages  of  the  HYPER 
are  clear — not  only  is  it  accurate  but  also  it  reduces  the  execution  time 
tremendously.  Thus,  the  HYPER  should  be  considered  to  be  a  practical  high 
frequency  model. 
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7.  A  VARIABLE  OENSITY  PARABOLIC  EQUATION 


Gregory  A.  Krelgsmann 

Northwestern  University 

Ding  Lee 

Naval  Underwater  Systems  Center 

Frederick  Tappert 
University  at  Miami 


ABSTRACT:  In  this  paper,  we  derive  a  new  parabolic  equation  (PE) 
that  Incorporates  the  effects  of  a  variable  ocean  density.  This 
density  can  be  smooth  or  plecewlse-smooth.  Thus,  our  model  reduces 
to  the  standard  PE  when  the  density  Is  constant  and  It  alleviates 
the  need  for  Interfaclal  conditions  when  the  density  Is  stratified 
In  a  piecewise  fashion.  We  also  present  a  numerical  scheme  that 
will  be  used  to  solve  our  new  equation.  This  difference  scheme  has 
a  conservation  law  that  Is  the  discrete  analog  of  the  new  PE's 
conservation  law. 
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DERIVATION 

The  propagation  of  sound  in  an  ocean  with  variable  density  p  is  governed 
by  the  elliptic  equation 

pv  •  j  P  +  k'2n2p  »  0  ,  (1) 

where  p  is  the  acoustic  pressure,  k*  «  w/cQ,  u  is  the  frequency  of  the  time 
harmonic  source,  cQ  is  a  reference  sound  speed,  n  -  c0/c,  ano  c  is  the 
sound  speed  in  the  ocean.  A  time  dependence  of  e“ilut  is  suppressed. 

Equation  (1)  is  to  be  solved  in  a  spatial  domain  D',  which  contains  the 
water.  A  simple  model  is  obtained  by  assuming  both  the  ocean  bottom  ano  the 
water-air  interface  are  flat.  Specifically, 

D'  -  (x'.y'.z*)  |  x '  |  <- ,  |y'|  <  «,  0<z‘  <H' 

where  the  primeo  variables  denote  dimensional  quantities.  Since  equation  (1) 
is  elliptic,  boundary  conditions  are  required  to  complete  the  mathematical 
description  of  the  problem.  The  conditions  used  in  this  report  are 

f|r  -  0,  z1  »  H1  (2) 

and 

p  =  0,  z '  *  0  (3) 

Thus,  the  ocean  has  a  hard  bottom  ana  a  pressure  release  (free)  surface. 
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The  source  deriving  equation  (1)  is  usually  modeled  as  a  point 
disturbance  located  at  x‘  *  y1  *  0,  z'  =  zQ'.  it  is  omitted  from  equation 
(1)  for  simplicity. 

In  many  underwater  applications  the  domain  (in  polar  coordinates),  O'  = 
|(r',z',e)  1 0  <_  r '  £  R ' ,  0  <  z'  H* ,  0  £  e  £  2wj  ,  where  equation  (1 )  must  be 
solved  is  extremely  slenaer.  By  this  we  mean  the  parameter 

e  -  (H'/R')2  (4) 

satisfies  e  <<  1  where  R'  is  the  maximum  range  of  interest.  We  now  introduce 
the  dimensionless  variables  r  and  z  used  by  Tappert  [1],  i.e.. 


r  -  ck'r'  (5) 

and 

z  «  k'z'  (6) 


Accordingly  D'  is  transformed  into 


D  = 


0  _<  r  <  z. 


0  <  z  0  l. 


0  <  ©  <  2n 


(7) 


where  a  =  (k'H'jH'/R*  .  We  asume  that  this  number  is  fixed  and  is  order  one 
with  respect  to  the  parameter  e.  Introaucing  this  change  of  variables  into 
equations  (1)  -  (3)  we  find  the  acoustic  pressure  satisfies 


2 

C 


r  .  i 

1 

r  pz  1 

LPrr  +  F  pr 

- D„  P„ 

p  r  rrJ 

+  £[pzz  r  pzJ 

(8) 


© 


© 
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ana 

|f  =  0,  z  -  i  .  (10) 

In  aaaition  to  the  boundary  conditions  (9)  ana  (10),  we  demano  that  p  is 
boundea  as  r  ■+  i . 

We  now  make  the  assumption  that  n^  deviates  slightly  from  a  constant 
and  takes  the  functional  form 

n2(x\y',z')  .  1  +  cf (r,z)  .  (11) 

The  constant  1  in  tnis  equation  is  arrived  at  by  taking  cQ  to  be  the  average 
of  c  throughout  D.  The  factor  e  in  (11)  demonstrates  the  weak  dependence  of  c 
on  depth  and  range.  (This  apparent  minor  perturbation  creates  profound 
effects  on  acoustic  propagation  when  the  range  is  as  short  as  a  few 
wavelengths i ) 

We  also  assume  that  the  density  p  depends  upon  the  variables  r  ana  z  in  a 
smooth  or  piecewise  smooth  fashion. 

When  (11)  is  insertea  into  (8)  we  observe  the  presence  of  the  small 
parameter  e  in  front  of  nearly  every  term.  To  cavalierly  set  these  terms  to 
zero  would  render  a  physically  meaningless  result.  Guided  by  previous 
experieme  with  such  matters,  we  apply  the  method  of  multiple  scales  to  this 
equation.  Specifically,  we  assume  that 

PU'.y'.z'  )  =  P(S,r,z;  e)  ,  (12) 
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where  the  fast  van  able  5  is  defined  by 

5  -  r/e  . 


(13) 


Inserting  this  variable  and  (11)  into  (8),  we  obtain  the  equation 

[p5{*p]  *  ‘K5  ~  pi  *  fzz  ~  ~7  pz  *  fp] 

*  '2[prr  *?pr  -TPr]  •  <‘«> 

The  subscripts  denote  partial  differentiation.  Next  we  make  the  assumption 
that  P  has  the  asymptotic  expansion 

00 

p  Pn(5,r,2,e)  ,  e-vO  .  (15) 

n=Q 

When  this  expression  is  inserted  into  (16)  we  equate  to  zero  the  coefficients 
of  the  powers  of  e.  This  yields  an  infinite  sequence  of  equations  of  which 
the  first  two  are 


L  P0  8  pOCC  +  p0  -  0 


(16) 


ana 


L  P, 


1 

?p  +  —  P _ L  p 

Ore  r  r0?  p  US 


+  P 


Ozz 


-  —  P  + 
p  oz 


fp. 


(17) 
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for  0<z<£,  0  <  r  <  a.  Inserting  (15)  into  the  boundary  conditions  (9)  and 
(10)  ana  equatinc  to  zero  the  coefficients  of  the  powers  of  e,  we  obtain  an 
infinite  sequence  of  boundary  conditions.  The  companions  for  (16)  and  (17)  are 

P„  «  0.  z  -  0,  n  -  0,1  ,  (18) 

and 

3P 

j~  «  0,  z  «  l;  n  .  0,1  .  (19) 

The  solution  of  equation  (16)  is 

P0  -  A0(r,z)  eU  +  BQ(r,z)  e-is  ,  (20) 

where  the  amplitudes  Aq  and  Bq  are  functions  of  tne  listed  variables. 

Because  of  the  assumed  time  dependence,  e~ita)t,  we  set 

B0(r,z)  -  0  (21) 

as  a  failure  to  do  so  would  yield  incoming  waves  from  infinity.  Inserting 
(2C)  ana  (21)  into  (17)  gives 


i_  P 


i  *  [: 


21  A0r  +  F  A0  "  p  A0  +  A0zz  -  T  A0Z  *  f  A«1  e<t  •  l22' 


'o] 


which  has  the  general  solution 


P:  -  Aj(r,z)  elC  +  j  M(Aq)  elF’ 


(23) 
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where  M(Aq)  is  the  bracketed  term  on  the  right-hand  side  of  equation  (22) 
We  observe  that  remains  bounded  as  5  ■  r/o«  only  if 

M(V  -  21  V  +  7  \  -  A0  +  A022  -  A02  +  f  A0  »  0  . 

If  we  now  set 


U 

Ag(r,z)  «  Vp(r,z)  ^ 


into  (24),  we  find  that  uQ  must  satisfy  our  new  variable  density  parabolic 
differential  equation  (VDPE) 


3  U 

-2i  *f  “0 


We  now  make  a  few  interesting  observations  about  this  new  equation. 

First,  when  p  is  a  constant,  (26)  reduces  to  the  standard  parabolic  equation 
[lj.  Secondly,  the  differential  operator  involving  z  is  symmetric  or  formally 
self-adjoint  [2].  Thirdly,  the  quantity 


,  /Kf* 


is  independent  of  range,  i.e.,  dE/dr  h  0.  This  follows  directly  fmm  (26)  ana 
the  boundary  conditions  for  uQ  are 


uo  =  0* 


0  , 
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ana 


fr^uo)-°-  z-‘  •  <29> 

Equations  (28)  and  (29)  are  direct  consequences  of  (18),  (19),  and  (25). 
Fourthly,  we  observe  that  equation  (26)  itself  was  derived  without  using  the 
bounaary  conditions  given  in  (18)  and  (19).  Thus,  our  new  parabolic  equation 
will  hold  even  when  more  realistic  boundary  conditions  are  implemented. 
Finally,  equation  (26)  can  be  used  even  when  p  is  piecewise  smooth.  This  will 
allow  us  to  study  interfaces  that  are  not  planar  or  straight  lines.  In  this 
sense  our  new  parabolic  equation  extends  the  analysis  given  by  Lee  and 
McDaniels  [3,4]. 


A  CONSERVATIVE  FINITE  DIFFERENCE  SCHEME 

In  this  section,  we  present  a  finite  difference  scheme,  which  is  second 
order  accurate  in  depth  and  first  order  accurate  in  range,  for  solving 
equation  (26).  This  difference  scheme  will  conserve  a  discrete  analog  of 
dE/ar  =  0,  where  E  is  given  by  (27).  The  method  of  analysis  and  other 
examples  are  given  by  Kriegsmann  and  Mahar  [5]. 

We  begin  by  rewriting  (26)  as 

-2i  TT  "  a  il  [b  JI  <auo>]  +  fuo  ’  <30> 
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where  d  ^  ana  b  h  1/P.  Setting  uj  =  u0(rn,Zj),  it  is  easily 
verified  by  Taylor's  theorem  that 


a  77  t15!!  (-0)] 


(rVo} 


L(uj)(az)“2  +  O(az)2  , 


(31) 


where  r  >  nar,  z.  ■  jaz,  and  L  is  defined  by 

I  i  J 


f/ ,,n\  an  .n  r  n  n  n  ,,ni  .  an  Kn  r  n  ,,n  n  ni 

L  uj  ■  a,,  b  1  a...  u..,  -  a,  u.  +  a.  b  ,  a.  ,  u.  .  -  a  u. 

J  j  j+1  I  J+1  J+1  j  Jj  J  j_l  L  3-1  j-l  J  Jj 


The  Crank-Nicholson  scheme  for  solving  (30)  is  the  one  we  shall  use.  It 
is  simply 


where  a  -  (1/2)  ar/(az)2,  0  -  ar/2,  and  Uj  is  the  numerical 
approximation  of  Uj.  Equation  (33)  is  solved  in  the  usual  fashion. 

We  now  define  the  vector  Un  by 


u" 


(34) 


where  the  superscript  T  denotes  a  transpose.  The  quantity  E  defined  by  the 
*>2  -  norm  of  U11,  i  .e. , 
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Q 


r.. 


^  : 


a 


lu 


n,  2 
i 1 


(35) 


is  the  discrete  analog  of  E  defined  by  (27).  We  shall  now  prove  that  E  is 
range  independent,  i.e., 


c 

fcn+l  "  En,  for  all  n. 
Defining  W,  by 

J 

wj  -  u3+1  +  “3  . 

ano  multiplying  (33)  by  Wj  we  obtain 


-2i 


UJ+1|  2  "  |  Uj|2  +  Rj  “  xHj  *LWj  *  0f"lW«  l2 


Jl  J 


The  term  Rh  rea]  and  gjven 


r"  s  -■«  k  °r‘]  • 


Summing  (38)  from  j  *  1  (zQ  *  o)  to  j  »  N  (zu  .  &  ),  we  obtain 


1-  ^  i  N  N 

• 

"2l  {En+1  En}  *  A  LWJ  +  8  j?0 

“j 

-  R. 


(36) 


(37) 


(38) 


(39) 


(40) 
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*  *>n 

where  Rn  is  the  norm  °f  the  vector  R,  with  components  R(j, 

defined  as  in  (34).  The  last  two  terms  on  the  right-hand  side  of  (40)  are 

N 

real.  The  result  given  in  (36)  follows  because  the  term  x  >.  B.  LW  ■  is  real 

>0  J  ^ 

also.  To  verify  this  fact  we  rewrite  this  sum  as 

N  N  2 

Z  "j  l“j  -  -So  *  Vi  -  Z  I  aj“j  -  ■  <«> 

j-0  j-0  J"7 


where 


V  ao  bl  Ho(al  W1  "  ao  Wo)  ’ 

7 


Vl  =  aN  bn+l  «n(Vi  WN+1  "  aN  Wn) 


Now  the  gQ  term  is  zero,  because  W0  »  Uq+*  +  u£  ana  both 

Uq+1  ana  u£  are  zero.  The  third  term  in  (41)  is  real.  The  term 

9n+I  vanishes.  This  is  because  the  bracketed  term  in  (43)  is  the  discrete 
implementation  of  the  boundary  condition  (29). 
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8.  CONCLUSIONS 

'Mathematical,  physical,  and  computational  contributions  to  the 
underwater  acoustic  wave  propagations  have  been  made  due  to  combined 
multiple  efforts  among  the  authors  presented  here.  These  accomplishments 
not  only  enhance  the  PE  capability  by  extending  the  solution  to  three- 
dimensional  problems,  but  the  work  can  be  extended  to  handle  acoustic  wave 
propagations  In  elastic  media.  It  Is  necessary  that  a  complete  numerical 
computation  be  performed  to  confirm  the  validity  of  these  theoretical  and 
computational  developments. 

^  loA  ojrrhi  -■>  A&  !  >  s  ' 
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